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The continuous improvement in localization errors (sky position and distance) in real time as LISA observes 
the gradual inspiral of a supermassive black hole (SMBH) binary can be of great help in identifying any prompt 
electromagnetic counterpart associated with the merger. We develop a new method, based on a Fourier decom- 
position of the time-dependent, LISA-modulated gravitational-wave signal, to study this intricate problem. The 
method is faster than standard Monte Carlo simulations by orders of magnitude. By surveying the parameter 
space of potential LISA sources, we find that counterparts to SMBH binary mergers with total mass M ~ 10 5 - 
10 7 Mq and redshifts z < 3 can be localized to within the field of view of astronomical instruments (~ deg 2 ) 
typically hours to weeks prior to coalescence. This will allow a triggered search for variable electromagnetic 
counterparts as the merger proceeds, as well as monitoring of the most energetic coalescence phase. A rich 
set of astrophysical and cosmological applications would emerge from the identification of electromagnetic 
counterparts to these gravitational- wave standard sirens. 
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I. INTRODUCTION 

One of the key objectives of the planned, low-frequency 
gravitational-wave (GW) detector LISA (Laser Interferomet- 
ric Space Antenna) is the detection of supermassive black hole 
(SMBH) binary mergers at cosmological distances. The ob- 
servation of these chirping GW sources would deepen our un- 
derstanding of (i) general relativity, e.g. by offering unique 
tests of spacetime physics in the vicinity of SMBHs [1-5], (ii) 
cosmology, by providing additional constraints on the lumi- 
nosity distance-redshift relation [6-8], (iii) large-scale struc- 
ture, by indirectly constraining hierarchical structure forma- 
tion scenarios [9-11, 13], and (iv) black hole astrophysics, 
e.g. by allowing accurate determinations of Eddington ratios, 
and other attributes of black hole accretion, in systems with 
SMBH mass and spin known independently, from the GW 
measurements [14—16]. 

From a purely astronomical point of view, one of the most 
attractive features of the LISA mission design is the possi- 
bility to constrain the 3-dimensional location (i.e. sky posi- 
tion and distance) of GW inspiral sources to within a small 
enough volume that the identification of potential electromag- 
netic (EM) counterparts to SMBH merger events can be con- 
templated seriously. Indeed, the accuracy of such LISA lo- 
calizations at merger are encouraging, with an error volume 
Sn x Sz = 0.3 deg 2 xO. 1 for SMBH masses m x = m 2 = 1O 6 M 
at z = 1, for instance [17]. In Ref. [15], we have shown that 
this accuracy may be sufficient to allow an unique identifica- 
tion of the bright quasar activity that may be associated with 
any such SMBH merger. 

Another possibility, examined here in detail, is to monitor 
the sky for EM counterparts in real time, as the SMBH inspiral 
proceeds. This is arguably one of the most efficient ways to 
identify reliably (prompt) EM counterparts to SMBH merger 
events, since the exact nature of such counterparts is a priori 
unknown. Using the GW inspiral signal accumulated up to 
some look-back time, ff, preceding the final coalescence, one 



already has a partial knowledge of where the source of GWs is 
located on the sky. Since the sky position is deduced primar- 
ily from the detector's motion around the Sun, one anticipates 
that angular positioning uncertainties will not change too dra- 
matically during the last few days before merger, so that a 
targeted EM observation of the final stages of inspiral may 
be a feasible task. Here, we present an in-depth study of the 
potential for such pre-merger localizations with LISA, while 
we discuss various astrophysical concepts and observational 
strategies for EM counterpart identifications in a companion 
work [18]. 

The main purpose of the present analysis is thus to deter- 
mine the accuracy of SMBH inspiral localizations with LISA, 
as a function of look-back time, ff , prior to merger. The LISA 
detector is not uniformly sensitive to sources with different 
sky positions and angular momentum orientations. Results 
will thus generally depend on the fiducial values of these an- 
gles. Our first objective is to calculate the time-dependence 
of distributions of localization errors, for randomly oriented 
sources, over a large range of values for the SMBH masses 
and source redshift. A second objective of our analysis is 
to estimate source parameter dependencies for these distribu- 
tions of localization errors, i.e. how the 3-dimensional (sky 
position and distance) localization error distributions depend 
on the fiducial sky position of GW sources. This is useful to 
understand which regions of the sky may be best suited for the 
identification of EM counterparts to SMBH merger events. To 
the best of our knowledge, this angle dependence has not been 
explored in detail before, not even in terms of final errors at 
ISCO (i.e. at ff = f; sco , when using the complete inspiral data- 
stream, up to the innermost stable circular orbit, or ISCO). 

Parameter estimation uncertainties for LISA inspirals have 
been considered previously, under a variety of approximations 
[5, 7, 8, 13, 17, 19-22]. These studies differ in the levels 
of approximation adopted for the GW waveform, using vari- 
ous orders of the post-Newtonian expansion. The LISA signal 
output for these waveforms are obtained through a linear com- 
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bination of the two GW polarizations, h+(t) and h x (t), with 
the beam pattern functions, F + and F x . The beam patterns 
define the detector sensitivity for the two polarizations. They 
are determined by the angles describing the instantaneous ori- 
entation of the LISA constellation relative to the GW polar- 
izations. As the LISA detector constellation orbits the Sun, 
with a one year period, F + and F x are slowly changing in time 
and this introduces an additional time dependence in the LISA 
signal. As first shown by Cutler [19], the source sky position 
can be determined with LISA using this modulation. In the 
formalism given by Cutler [19], this modulation couples time 
and angular dependencies in a complicated way, making the 
estimation of localization errors numerically costly for a large 
set of SMBH binary random orientations and parameters. 

Using a different approach, Cornish & Rubbo [23] have de- 
rived the orbital modulation in a much simpler form, in which 
the angular parameter dependence and the time dependence 
can be decoupled. Here, starting directly from the original 
Cutler [19] expression, we give an independent derivation of 
the Cornish & Rubbo [23] formula and write it in an equiv- 
alent form, from which decoupling is more evident. We do 
this by expanding the LISA response function into a discrete 
Fourier sum of harmonics of the fundamental frequency of 
LISA'S orbital motion, /© = lyr" 1 . Since LISA'S orbit does 
not include high frequency features, we expect this sum to be 
quickly convergent. In fact, it is clear from the Cornish & 
Rubbo [23] result that the expansion terminates at 4/© and 
that there are no higher order harmonics due to the detec- 
tor's motion. The series coefficients in the expansion are in- 
dependent of time and only depend on the relative angles at 
ISCO. We then develop a Fisher matrix formalism in which 
parameter error distributions can be mapped independently of 
time, while the time dependence can be computed indepen- 
dently of the specific SMBH binary orbital elements. A Monte 
Carlo simulation for random binary orientations then becomes 
a simple linear combination, without any integral evaluations. 
This greatly reduces the numerical cost of estimating param- 
eter uncertainty distributions, even at fixed observation time 
(e.g. to map distributions of errors at ISCO). We use this nu- 
merical cost advantage 

1 . to map the distribution of localization errors for the full 
three dimensional grid of SMBH total mass (M = 10 5 - 
1O 8 M ), redshift (z = 0.1-7) and arbitrary look-back 
time (ff) before merger, 

2. to study how source localization error distributions vary 
systematically with sky position, and 

3. to discuss implications, in terms of advance warning 
times, for prompt electromagnetic counterpart searches 
with large field-of-view astronomical instruments. 

We call this new approach the harmonic mode decomposi- 
tion (HMD). The method verifies that the amplitude modula- 
tion, which is restricted to frequencies less than 4/© = 1.3 x 
10" 7 Hz, is indeed a very slow modulation when compared 
to the GW frequency of LISA SMBH inspirals (0.03 mHz- 
1 Hz). One plausibly expects that physical parameters which 



determine the amplitude modulation (like the source sky po- 
sition and orbital inclination relative to the detector) can be 
estimated independently of the parameters which determine 
the GW frequency (like masses, orbital phase, time to ISCO). 
In the HMD method, the two sets of parameters are naturally 
separated and can be estimated independently. In particular, 
parameters related to the modulation can essentially be deter- 
mined on a background of GW-cycle averaged signal. In the 
present work, we compute LISA inspiral localization errors 
with the approximation that high frequency signal parameters 
have strictly no cross-correlations with parameters related to 
the slow orbital modulation. In addition to the numerical ad- 
vantages mentioned above, the HMD formalism offers a clear 
interpretation of the time evolution of uncertainties for the 
slow modulation parameters. This can be used to gain a better 
understanding of the general evolutionary properties of local- 
ization errors. The following questions, that we address in 
detail in our work, are particularly relevant. 

(i) Under what conditions do the localization uncertainties 
scale simply with the measured signal-to-noise ratio, 
and how do these uncertainties evolve during the final 
stages of inspiral? 

(ii) To what extent do the high and low frequency signal 
parameters decouple? 

(iii) What are the best determined combinations of the an- 
gular parameters? 

(iv) How and why does the shape of the 3D localization er- 
ror ellipsoid change during the final week(s) of obser- 
vation? 

In our analysis, we neglect the "Doppler phase" due to 
LISA'S orbital motion, SMBH spin precession effects and any 
finite SMBH binary orbital eccentricities. These approxima- 
tions are advantageous for the resulting simplicity, but the 
use of the HMD method is not restricted to these approxi- 
mations. We also outline a generalized HMD method which 
remains numerically much more efficient than standard meth- 
ods. We leave a numerical implementation of this general 
HMD method to future work. It will be particularly interest- 
ing to determine how our approximate results for the evolution 
of LISA localization errors are modified when spin precession 
effects are included, since spin precession effects were shown 
to improve the final localization errors by factors of 3-5 at 
ISCO [17, 22]. 

The remainder of this paper is organized as follows. In § II 
we define our conventions and the assumptions made in our 
analysis. In § III we expand the LISA GW signal in Fourier 
modes and obtain the conversion from actual physical parame- 
ters to corresponding Fourier amplitudes. In § IV, we incorpo- 
rate these results into a Fisher matrix formalism and derive the 
expressions necessary to estimate correlation errors for HMD 
signals. In § V, we quantify the computational advantages of 
the HMD method. In § VI, we present results from Monte 
Carlo computations of the time evolution of localization er- 
rors and discuss results in terms of advance warning times for 
prompt electromagnetic counterpart searches. In § VII, we 
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develop toy models to interpret the time-dependence of LISA 
localization errors and to answer questions (i)-(iv) above. We 
summarize our results and conclude in § VIII. 



II. ASSUMPTIONS AND CONVENTIONS 

This section is divided into three parts. First, we list the 
definitions of physical quantities used in this paper, in partic- 
ular the variables describing a SMBH inspiral. Second, we 
give the equations which determine the LISA inspiral signal. 
Third, we state all the assumptions made in this work. 



A. Definitions 

In general, an SMBH inspiral is described by a total of 
17 parameters. These include 2 redshifted mass parameters, 
(A4 z ,r] z ), 6 parameters related to the BH spin vectors, p sp m, 
the orbital eccentricity, e, the source luminosity distance, d^, 
2 angles locating the source in the sky, (On, fa), 2 angles that 
describe the relative orientation of the binary orbit, (8nl, 4>nl), 
a reference time, f mer ger, and a reference phase at ISCO, faco, 
and the orbital phase, fat,. Throughout this work, we restrict 
ourselves to circular orbits by omitting the orbital eccentricity, 
e, and instead of the orbital phase, fa b , we use the look-back 
time before merger, t , as our evolutionary time parameter. The 
LISA signal for a GW inspiral is determined by the above set 
of parameters and two additional angular parameters describ- 
ing the orientation of LISA, (S, $). We elaborate on the defi- 
nitions of our mass and angular parameters below. 



1. Mass Parameters 



merger is (e.g. eq. 3.3 inref. [25]) 

s 3/8 

MM z ,t)=*—r^Mf\ 

= 2.7xlO- 4 Hz^j % 3 /> 6 f 8 , (1) 

or equivalently 

t (M z ,f) = 5(8nf)-^ 3 Mf 3 

/ f \ "8/3 

= 6.7rnin(^J %25^& > ( 2 ) 

where Mg z is the redshifted total mass in units of 4 x 10 6 M Q , 
770.25 = 77/O.25 is the symmetric mass ratio (770.25 = 1 for equal 
component masses, § IIC),/ c = c/R® = c/(l AU) = 2.00mHz 
is the inverse light-travel time across the radius of the LISA 
orbit, and the null index stands for the order of approximation. 
The inspiral phase extends until the innermost stable circular 
orbit (ISCO), at 6M, is reached 

/ < /isco = 6" 3 / V'M; 1 = 1 . 1 mHz x Af£ , (3) 
f>fisco = 5(3/2)V 1 ^ = 33minxr/o 1 25 ^. (4) 

where fisco is the (observer-frame) look-back time before 
merger corresponding to the ISCO, and /isco is the (observer- 
frame) frequency at ISCO. 

In the present work, we fix the start of the observation 
(i.e. when the source first enters LISA'S frequency band) at 
look-back time f;, and examine how the value of an end-of- 
observation time, tf, prior to merger affects the precision on 
source localization. We restrict ourselves to pre-ISCO inspi- 
ral signals, corresponding to tf > fisco- Note that any instanta- 
neous look-back time t associated with an observation lasting 
from look-back times t\ to tf must obey fisco < tf < t < ti in 
our notation. 



For component masses mi and rm, the total mass is M = 
mi + m.2, the reduced mass is /x = m\tnilM, the symmet- 
ric mass ratio is 77 = \xjM and the chirp mass is defined as 
M. = Mifl 5 [24]. Throughout this work, we use geometrical 
units: G = c = 1 . In this case, the mass can be expressed in 
units of time: 10 6 M Q = 4.95 sec. The measured GW wave- 
forms are insensitive to the cosmological parameters, if they 
are expressed in terms of the luminosity distance and the red- 
shifted mass parameters, e.g. m z = (l+z)m (same for red- 
shifted chirp and reduced masses). 



2. Time Parameters 

We write a generic look-back time (or "observation time") 
before merger as t , and a generic redshifted GW frequency (or 
"observation frequency") as /[37]. We use the leading order 
(i.e. Newtonian) approximation for the frequency evolution. 
Therefore, the observed frequency at look-back time t before 



3. Angular Parameters 

LISA is an equilateral triangle-shaped interferometer with 
an arm-length of 5 x 10 6 km, orbiting around the Sun. The 
constellation trails 20° behind the Earth and is tilted 60° rel- 
ative to the ecliptic. The detector plane precesses around the 
orbital axis with the same one-year period as the orbital period 
[26]. 

Following closely Refs. [19] and [17], including in nota- 
tion, we define two coordinate systems. The barycentric frame 
is tied to the ecliptic, with x,y lying in the ecliptic plane and 
z normal to it. The detector reference frame tied to the de- 
tector, with z' normal to the detector plane, while x',y' are in 
the plane and co-rotating with the detector so that the arms 
are described by time-independent vectors. We refer to the 
barycentric frame with normal coordinates and to the detector 
frame with primed coordinates. The unit vectors defining the 
source location on the sky, N, and the SMBH binary orbital 
angular momentum, L, are described by polar angles (0 N , fa) 
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and {6 L ,<f>L) in the ecliptic frame, (6' N ,(p' N ) and (9' L ,4>' L ) in the 
detector frame: 

N(0n, 4>n) = zcos#Ar + xsin6>AiCos0Ai + ysin6>;vsin</>;v, (5) 
L(^l,^l) = zcos# L + xsin# L cos0z, + ysin# L sin0 L . (6) 

Since we assume no SMBH spins, orbital angular momen- 
tum is conserved and the (0 N , <j> N , L , <j> L ) coordinates are time- 
independent properties of the sources. 

Let (S, $) be the two angles specifying the orientation of 
the LISA system in the ecliptic: $ describes its orbital phase 
during its motion around the Sun, while E describes the ro- 
tation of the triangle around its geometrical center. If their 
values at merger are written So and <I>o, then at an arbitrary 
look-back time f: 



E(t) = Eo-(Js,t, 

$(f) = $ -Wffif, 



(7) 
(8) 



where = 2tt/ yr is the orbital angular velocity around the 
Sun. 

The time dependence of the detector normal vector z' can 
be expressed as 



z = -z xcos<P ysin<P. 

2 2 2 J 

The detector angles are given by 

1 V3 
cos6' N = -cos6 N — — sin^cos($-^>jv), 



//„ = £+tan 1 



^ cos 6 N +\ sin 6 N sin($ - <f> N ) 
sin On sin($ - fa) 



(9) 

(10) 
(11) 



Let us also define tjj', the polarization angle of the GW wave- 
form, as [17] 



tan?// = 



L-z'-(L-N)(z'-N) 
N • (L x z') ' 



(12) 



Note that there are only 6 independent angular parameters 
(6nf,<f>N,9L,<pL,E,$). Other detector specific quantities like 
6' N , <fi' N , 9' L , <p' L , and ip' can be expressed in terms of these 6 
independent parameters using eqs. (5-12). 

Let us introduce a new set of 6 independent angles, 



H = (0jv, 4>n,9nl, <(>nl, a, 7)j 
with the following definitions: 



(13) 



• 6>;vl is the relative latitude of L and N (i.e. the inclina- 
tion of the binary orbit to the line of sight), 

• (j>NL is the relative longitude of L and N, 

• a = E-^ + fa-^f, 

• j(t) = $(t)-fa. 



The explicit definitions are given in Appendix B. 

Let us refer to the angles at the reference time t = f mer ger = 
as £1(0). Although $ = $(f) and E = E(t) are time-dependent, 
as given by (7,8) a is a time-independent combination, unlike 
the time-dependent 7 = 7(f). The angles at t = are thus given 
by 9.{0) = {6 N ,4> N ,e NL ,4> NL ,a,^ ). 

These angles have the interesting property that they pos- 
sess isotropic a priori distributions, like the original i!(0) vari- 
ables, but the measured GW waveforms expressed in terms of 
these new variables are much simpler than when they are ex- 
pressed in terms of the original set eqs. (5-12). 

Two additional quantities which are useful to describe the 
sensitivity of the detector in various directions are the antenna 
beam patterns [19]: 



Px,+' 



± cos 6' N sin 2(j)' N sin 2ip' N , 



(14) 



where the sign ± is defined to be positive for F x , and negative 
for F+. Equation (14) and the transformation rules eqs. (5-12) 
define the time evolution of the antenna beam patterns for a 
given set of final angles £1(0) as the LISA system orbits around 
the Sun. Note that the LISA system is equivalent to two inde- 
pendent orthogonal-arm interferometers which are rotated by 
45° relative to each other [19]. Both data-streams are given 
by the same equations (see eq. [21] below), modulo a change 
of one of the angles for the second detector: <j)' N = <p' N -ir/4 
(or equivalently a ll = a 1 -n/4 using our time-independent an- 
gular variables). Thanks to this simple relationship between 
the two data-streams, it is possible to carry out all the calcu- 
lations for the first data-stream, and later include the second 
data-stream in the final expression by varying the fiducial an- 
gle a. 



4. Grouping the Parameters 

We group the most important parameters describing the in- 
spiral as follows: 



Pslow = {<4,^}, 



Wast 
Pspin 



(15) 
(16) 

{2 spin magnitudes, 4 spin angles}. (17) 



{M z j Mz 7 ' merger j 0isco}, 



This organization of parameters has fundamental importance 
in our formalism. As we show in § II B, the parameters pf ast 
and p S pi n relate to the high frequency GW signal, while the pa- 
rameters p s i ow relate to the distinctly slow orbital modulation. 



B. LISA Inspiral Signal Waveform 

For a circular binary inspiral, the two polarizations of GW 
signal are well approximated by the restricted post-Newtonian 
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expressions 



h+(t) = 2 , JJ (l + cos 2 gjv L )cos(few(0, (18) 



MO 



d L 

^ M 5/3 (7Tf} 2/3 



cos 6nl sin 0gw(O ■ (19) 



The GW phase </>gw(0 = <AGw(f>fast,Pspin;0> which is twice the 
orbital phase , <p(t) = 2(f> OI \,(t), can be expanded into the series 

0Gwfeast,Pspin;O ~0isco + (fo(M z ; f ) + </>i(M z , fj, z ; t) 

+ MM z ,fj, z ,p S p a ;t) + ..., (20) 

where <fio(M z ; t) is the leading order Newtonian solution to the 
phase evolution, successive terms correspond to small general 
relativistic corrections, (frsco is the reference phase at ISCO 
and <p„(tisco) = for all n > 0. The instantaneous GW fre- 
quency is defined as the time derivative of the GW phase (20), 
i.e. / = f(t) = dcf>cw /At, which changes very slowly compared 
to the GW phase itself, (/>gw(0- In practice we use the New- 
tonian approximation (1), fo(t) = dfo/dt. Note that equation 
(20) depends implicitly on the reference time, f mer ger> since 
our time variable t is the look-back time before f mer ger (see 
§HA2) 

The signal measured by LISA is a linear combination of the 
two polarizations (18), weighted by the antenna beam patterns 
F+ u and f' ,h for each of the two equivalent interferometers, 
defined by (14), resulting in the two observable data-streams 



h^(t)=^[F l /%(t) + F^% 



(OL 



(21) 



where the factor \/3/2 = sin(60°) comes from the opening an- 
gle of the LISA arms. The beam patterns are determined by 
the relative orientation of the source polarizations and the de- 
tector. Their time-dependence is due to the following three 
main effects: LISA changes its orientation as it orbits the Sun, 
LISA changes its relative distance to the source as it orbits the 
Sun, and the orbital plane of the SMBH binary can precess be- 
cause of spin-orbit coupling effects. Substituting (18) in (21) 
and expressing it in complex form, we get 



MM Z J) G i, n ^Q ^ yy 0cw(pra St , Psp i„ ;0 ) (22) 
d L 



where A(M z ,f)/d^ defines the overall amplitude scale, with 

A(M z ,f) = 2V3(irf) 2 / 3 M 5 z /3 . (23) 

The G(Q,/) factor defines the angular dependence of the sig- 
nal, 



(24) 



where Ga(O), the amplitude modulation, captures the vary- 
ing detector sensitivity with direction and polarizations of the 
GWs, 

1 + c °S 2 #iVL r I.n /n> • a rUI/riN /OC n 

G A (fi)= 2 F+ (^-'cos^^x ( 25 ) 



The additional tpv(Cl,f) modulation is the Doppler phase 
modulation, which is the difference between the phase of the 
wavefront at the detector and at the bary center [19]: 



f 

(p D (fl,f) = 2tt— sin#ivCOS7. 

Tc 



(26) 



There is a non-negligible number of Doppler phase cycles 
only for a GW frequency satisfying / > f c (see definition 
of f c below eq. [2] above). However, equation (3) shows 
that / < /isco < fc, hence the f c frequency is reached only 
after ISCO for typical SMBH component masses of mi = 
ni2 = 1O 6 M and redshift z = 1. Even for smaller 10 5 M Q 
component masses, the total number of cycles, A^ pm , remains 
< 1 until the final 5hr of inspiral. Therefore the Doppler 
phase (26) is practically negligible for SMBH inspirals. In 
fact, estimating localization errors without accounting for the 
Doppler phase affects results by less than a factor of 10~ 3 (for 
ni\=ni2 = 10 6 M Q at z = 1; S. A. Hughes, private communica- 
tion). Therefore, in eq. (24), we neglect tp D (Q,f) and restrict 
ourselves to the approximation 



,i,n. 



(27) 



The explicit frequency-dependence dropped out, and the time 
evolution of the signal Ga is now fully determined by the time 
evolution of the angles fL 

Note that the amplitude modulation (25), G^ n (fi), is tra- 
ditionally expressed in complex polar notation (e.g. [19]), 
where the magnitude and argument of the complex number 
are called polarization amplitude and phase. As we will show, 
the mode decomposition is simplest in the original Cartesian 
complex form (25), which already includes both the polariza- 
tion amplitude and phase; thus, we do not distinguish these 
two quantities in the following. The function G 1,ll (fl,f) given 
in (24) also accounts for spin-orbit precession if the orbital 
orientation (9nl, <Pnl) in ^ is chosen to be time-dependent, to 
satisfy the equations for spin-orbit precession, and if an extra 
precession phase shift, e\p(iS P (6 NL , 4>nl)), is introduced (see 
eq. 2.14 in Lang & Hughes [22]) in addition to the Doppler 
phase in (24). In our calculations, we neglect spin precession 
but discuss how the HMD method can be extended to include 
that effect in§IVC. 

Finally, we express the measured signal (22) as 



h hll (p;t) = /z c (pfast,p S pin;f) x /& u (p s iow;0 



1.11/ 



(28) 



where h c is the high frequency carrier signal and h m is the 
slow modulation: 



fcc(Pfast,Pspin;0 = ^./(O)^'^*^" (29) 

«m (Psi ow ;0 = - J2 — 1 • (30) 



Equation (28) shows that the two sets of parameters p s \ ow 
and {pfasbPspin} are exclusively determined by the low fre- 
quency modulation and the high frequency carrier, respec- 
tively. For this reason, we only expect a low level of cross- 
correlation between these sets of parameters: parameters asso- 
ciated with very different timescale components should essen- 
tially decouple. In Sec. VII A and Appendix A, we consider 
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several toy models which allow us to understand the neces- 
sary conditions, and the extent to which, parameters associ- 
ated with high and low frequency components decorrelate in 
the course of an extended, continuous observation. 

C. Simplifying Assumptions 

In the present work, we make the following assumptions: 

1 . We assume that the amplitude modulation can be used 
to determine the luminosity distance and angular pa- 
rameters, psiow = {d L ,9 N , 4>n,&nl, 4>nl}, while the other 
parameters, p fast = {M z , fi z ,t meiger ,(j)isco}, are deter- 
mined from the high frequency GW phase. We assume 
no cross-correlations between these two sets of param- 
eters. This is supported by the results listed in Table 1 
of Hughes [7], which shows the full co variance matrix 
of a Monte Carlo realization of 2PN waveforms. The 
correlation coefficients are ^0.1 for the above quanti- 
ties, and the absolute scale of the second set of param- 
eters is very low in the first place. Berti et al. [12, 13] 
also report that the sets p fast and p slow are relatively un- 
corrected for general relativity and even for alternative 
theories of gravity. In the latter case, the carrier h c {t) in 
the signal (28) is modified but not the slow modulation, 
h m (t), so that the general expectation of decoupling is 
maintained. 

2. We assume that there are no additional errors on the 
detector orientations $(0) and 5(0). These parame- 
ters are given by f merg er via eq. (8) and (7), and / me r ge r 
itself is determined by the high frequency carrier sig- 
nal to high precision. Using the full data-stream up to 
ISCO, <5f merg er ~ 2 sec [5, 7]. Using (8) and (7), we es- 
timate | <5$(0)| = 1 55(0) | = w©<5f me r g er = 4 x 10" 7 rad = 
0.08". This is so small that we expect the errors 
£$(0) and (55(0) to be negligible at any relevant end-of- 
observation times tf > fisco, even if the ff-dependence 
of these errors scale as steeply as (S/N)' 1 (see also Ap- 
pendix A). 

3. We use the circular, restricted post-Newtonian (PN) ap- 
proximation for the GW waveform, keeping only the 
leading order (i.e. Newtonian) term in the signal am- 
plitude. Higher order corrections to the GW amplitude 
introduce additional structure to the waveform. They 
improve the parameter estimation uncertainties for high 
mass binaries [27, 28] and introduce additional corre- 
lations between the parameters. It will be important to 
consider these corrections to the amplitude in future in- 
vestigations. Arbitrary PN corrections to the GW phase 
only enter via h c in the signal given by eq. (28). Since 
we neglect correlations between the sets of parameters 
p s i ow and (pfast,Pspin), all the restricted PN corrections to 
the phase drop out and become irrelevant for the p s i ow 
parameter estimations. 

4. We neglect the effects of Doppler phase modulation. 
This is plausible for SMBH binaries with component 



masses m > 1O 5 M , since the GW wavelength in this 
case is generally greater than LISA'S orbital diameter 
and N pm < 1 (see eqs. [26] and [1]). 

5. We neglect SMBH spins and, in particular, neglect the 
spin-orbit precession for angular determinations. This 
assumption is useful in simplifying our equations and 
in focusing on the behavior of pure angular modulation. 
Future studies can incorporate spin-orbit precession by 
convolving the angular modulation decomposition with 
spin-orbit effects. 

6. We fix the start of LISA observations at a look-back 
time f; = min{fo(/mi n ), 1 yr} prior to merger. This cor- 
responds to the time when the GW inspiral frequency / 
crosses the low frequency noise wall of the detector at 
/min = 0.03 mHz, but we limit the initial look-back time 
to a maximum of 1 yr before merger. Note that LISA'S 
effective mission lifetime is estimated to be 3 yr. Inte- 
grated observation times longer (but also shorter) than 
our assumed t\ values are possible in principle, depend- 
ing on source specifics. In a more elaborate treatment, 
one could define t\ as an a priori random variable. We 
fix fi here mostly for simplicity and focus on the effects 
of varying the values of tf (< f;). In § VII A we show 
that localization errors are primarily determined by the 
end-of-observation time, tf, and that values of f; > 1 yr 
do not significantly change the evolution or final local- 
ization error estimates. If, however, t\ <C 1 yr (that is, 
if f merge r is within a few months of the beginning of ob- 
servation), then localization errors can become signifi- 
cantly worse than in our results with t\ = 1 yr. 

7. We neglect finite arm-length effects and we do not make 
use of the three independent observables of the time de- 
lay interferometry [29]. This is a valid assumption for 
SMBH inspirals since here / < c/L = 0.01 Hz. 

8. We neglect any finite orbital eccentricities. We note 
that, for eccentric orbits, higher order harmonics ap- 
pear in the GW phase. In principle, since these har- 
monics affect the high frequency GW phase, but not the 
slow modulation, including finite eccentricities should 
not significantly affect localization error estimates. For 
rather eccentric orbits, high-order harmonics with / ^> 
f c can have a non-negligible Doppler phase (2), which 
would lead to an improvement in the determination of 
9n and fa. Although eccentricity is efficiently damped 
by gravitational radiation reaction [30], the presence 
of gaseous circumbinary disks could lead to non-zero 
eccentricities for at least some LISA inspiral events 
[31,32]. 

9. We follow Barack & Cutler [21] in calculating the LISA 
root spectral noise density, S n (f), which includes the in- 
strumental noise as well as galactic/extra-galactic back- 
grounds. For the instrumental noise [13], we use the ef- 
fective non-angularly averaged online LISA Sensitivity 
Curve Generator[38], while we use the isotropic formu- 
lae for the galactic and extra-galactic background [21]. 
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10. Our analysis focuses on statistical errors and does not 
account for possible systematic errors. For example, 
waveform templates might be inaccurate either due to 
the imprecision of the theory if the true waveform is not 
the one predicted by general relativity, or due to practi- 
cal limitations from necessarily finite realizations of the 
large template space. Such inaccuracies can introduce 
new systematic errors. 



in. HARMONIC MODE DECOMPOSITION 

In our formalism, the angular information of the LISA in- 
spiral signal is contained exclusively in the periodic modula- 
tion due to the detector motion around the Sun, which adds an 
amplitude modulation to the high frequency waveform. This 
modulation has a fundamental frequency, /© = 1 / yr, along 
with upper harmonics jf^, where j is an integer. Although it 
is intuitively clear that the high frequency harmonics will tend 
to have a vanishing contribution, it is hard to establish this just 
by looking at eqs. (5-12), which define the time evolution. In 
this section we show that it is possible to derive surprisingly 
simple analytical expressions for the amplitude of each har- 
monic. We provide an outline of the derivation starting from 
the commonly used Cutler [19] formulae (5-12) and alterna- 
tively from those in Cornish & Rubbo [23]. We show that the 
derivation is much simpler in the latter case, in the sense that 
the Cornish & Rubbo [23] expression is almost already in the 
desired form. 



A. Derivation using Cutler [19] 



We expand the modulating signal (25,30) in a Fourier series 



MPsiow(0);f) : 



ddz) 



J2 s>(Psiow(0)y 7 ^', (3i) 



j=-oo 



where g 7 (p s iow(0)) are the mode amplitude coefficients and 
Psiow(O) are the distance and angle variables at t = (see 
§ II A). The coefficients can be obtained as 



£/(Psiow(0)) 



lyr 



d/G A (^(0)e" ,>ffl '. (32) 



Substituting the definition of G\(ft(t)) from eq. (25), using 
the time evolution of Q(t), eqs. (5-12), integral (32) can be 
evaluated. 

Although conceptually simple, a direct analytical evalua- 
tion of integral (32) is overly cumbersome. Thus, for practical 
reasons, we follow an alternative path. We start with the orig- 
inal Cutler [19] formulae, given by eqs. (14) and (25). First, 
using general trigonometric identities, we can express cos2x 
and sin 2x with tan(x) for x = (f>' N and x = ip' . In the second step, 
we express and substitute for tan^ and tan^ with ecliptic 
variables using (11) and (12). In the third step, we express the 
trigonometric functions in complex form. After this step, each 



term in the beam pattern (14) is of the form 



(33) 



where the sums over n and m integers are finite, containing 
only a few terms, and a„ and b„ depend only on the angles 
(On, 4>nl, ct). In the fourth step we simplify the product of frac- 
tions. It turns out that, after combining terms, the denomina- 
tors drop out exactly, leaving a formula just like (31), except 
that the largest element in the sum is \j\ = 4. In the fifth step, 
we substitute in (25), and finally, change back from complex 
to trigonometric notation for the coefficients, using the half- 
angles 0^/2 and 6 NL /2. Finally we arrive to the remarkably 
simple form: 

4 

MPsiow) = W 1 Y,(LNjDj + L*N*Dj + L*NjD* +LN*D*), 
j=o 

(34) 

where the functions L, N, and D depend only on the angular 
momentum angles, sky position angles, and detector angles, 
respectively: 



L(0, L ,>, L ) = sin'f^-'" 



/ / On \ . 4 _; / U N 

Dj(a,j) ee ie 2ia e ih 



(35) 



Nj(0 N ) ee WjCOS > [ I sin*"' I y ], (36) 



(37) 



where Wj = 1/16 x (9, 12\/3, 18,4^3, 1) for ; = (0, 1,2,3,4), 
respectively, and we have defined asterisks to refer to the fol- 
lowing conjugates: 



L*(0nl,4>nl) = L(tt-8 nl ,-4>nl), 
N*(0 N ) ee (-iy Nj (ir-9 N ), 



(38) 
(39) 



D*(a, 7 ) ee -D/-a,-7) = D/a,7). (40) 

Note that using these conjugate functions, only the non- 
negative terms < j < 4 remain in the sum (34). 

Substituting the time dependence implicit in 7 ee 7(0)+w©f , 
equation (34) becomes 



fcm(Pslow(0),0 = ^L(Z)" 1 ^-(PslowWy 7 "®', (41) 

where the coefficients are 

_ f LNjDy^ + ^NJDjiO) if j > 
Sj(Psiow(0)) - I ^N^D^iOy + LN^D^iO) if j < (42) 

and the detector functions Dj(0) and D*j(0) refer to their 
values at t = 0, 7(0). (Note that L,Nj,L*,Nj are all time- 
independent.) Since the decomposition (31) is unique, the 
coefficients (42) that we read off from our result also satisfy 
eq. (32). 
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B. Derivation using Cornish & Rubbo [23] 

Our result in (34) can also be derived from the Cornish & 
Rubbo [23] formulae for the LISA response function. In their 
paper, these authors use a different set of angles, which relate 
to ours as follows: 9 CR = 6 N , <p CR = fa, -0 CS = fa L -(n/2), 
\ CR = -a+fa = $-S-(37r/4), and a CR = -y+fa = $. Note 
that our set of angles is very similar to theirs, except that we 
measure the detector angles relative to the source, fa. This is 
advantageous given the rotational symmetry around the Earth 
orbital axis, making angles relative to the source the only ones 
that should have an effect on the measured GWs; we expect 
fa to drop out of the equations when using a and 7. Note, 
once again, that the variables (0 Nl fa,8 NL} fa L , a) are time in- 
dependent, while 7 = 7(f). Writing the Cornish & Rubbo [23] 
beam patterns for low frequencies, which is fully equivalent 
to Cutler [19], with our angular variables [39], we get 

= -I [cos(2fa L )D l f(t) -sin(2fa L )Df(t)], (43) 
F 1 / = -^[ & m(2fa L )D 1 /\t) + cos(2fa L )Df(t)], (44) 



where 



3^ {-36 sin 2 6 N sin(27+ 2a) + (3 + cos2^) 

x {cos(2</>at) [9 sin(-2a + 2 fa) - sin(47 +2a + 2 fa)] 

+ sin(2(/> A ,)[cos(47 +2a-2fa)-9 cos(-2a + 2 fa)] } 

- 4V3sin(26' iv )[sin(37 + 2a)-3sin(7 + 2a)]}, (45) 



and 



^ {\/3cosfl A r[9cos(-2aO-cos(47 + 4aO] 
- 6sin6W[cos(37+2a) + 3cos(7 + 2a)]}. (46) 

One notices instantly that the time dependence here is much 
simpler than in the original Cutler [19] formula, as it is in- 
scribed only in the various harmonics of 7. We can identify 
the highest harmonic present to be 47. Expanding the trigono- 
metric functions using standard identities, we obtain 



D + = —_ 



32 



/ 9(3 + cos26W) \ 
-12\/3sin2#Ar 

36 sin 2 On 
4V3sin26> w 
\ 3 + cos26W / 



and 



D x =-- 



/ -9cos6W \ 
6\/3sin0Ar 


2\/3sinftv 
cosftv 



I sin(2a) \ 
sin(2a + 7) 
sin(2a + 27) 
sin(2a + 37) 

V sin(2a + 47) / 



I cos(2a) \ 
cos(2a + 7) 
cos(2a + 27) 
cos(2a + 37) 

V cos(2a + 47) ) 



(47) 



(48) 



where a • b = J2 n a nbn is the usual dot product. Now, the sec- 
ond sets of elements carry the time dependence and the de- 
tector orientation information, while the first sets describe the 



sky position. Note that the explicit fa dependence dropped 
out, as expected. Next, we manipulate equations (47,48), sub- 
stituting complex expressions for the trigonometric ones, and 
substituting these into eq. (25). We finally arrive at eq. (34) 
after changing to half-angles 9 N /2 and 9 NL /2. 

We note that eqs. (34) or (41,42) are fully general ex- 
pressions, equivalent to the standard LISA inspiral signal in 
eqns. (14) and (25). The two data-streams are obtained by 
substituting a = a 1 and a", corresponding to the two indepen- 
dent LISA-equivalent Michelson interferometers (see § II A). 
To verify our final result, we compare numerically the signals 
computed using eqs. (14,25) with the signals computed us- 
ing eqs. (41,42), for a large set of random choices of angles. 
Agreement is achieved at machine precision levels. 

The main utility of eq. (34), is that it can be used to "de- 
construct" parameter error histograms, i.e. to understand how 
the errors depend on the fiducial values of the parameters. As 
compared to Cornish & Rubbo [23], our result leads to an ex- 
plicit decoupling of the signal angular dependence into simple 
products of one-dimensional functions. In particular, the de- 
pendence on sky position, angular momentum, and detector 
angles are separated. Using the special conjugate functions 
L* , D* , N* , eq. (34) displays the symmetry properties of the 
signal. Finally, one angular variable, fa is eliminated exactly. 



IV. ESTIMATING PARAMETER UNCERTAINTIES IN 
THE HMD FORMALISM 

Parameter estimations for LISA GW inspiral signals are 
possible with matched filtering and the expected uncertainties 
can be forecast using the Fisher matrix formalism [33, 34]. In 
this section, we apply this approach to the LISA signal derived 
in § III, with an angular dependence of the signal decomposed 
into harmonic modes. In § IV B, we consider the simple case 
of a high frequency carrier signal that is modulated by a low- 
frequency function, without any cross-correlation between the 
two sets of relevant parameters. We derive a simple formula 
for the estimation of modulating parameter uncertainties. In 
§ IV C, we consider a more general post-Newtonian signal and 
show that parameters related to source localization can still be 
decoupled from the time evolution and the other source pa- 
rameters. 



A. Fisher Matrix Formalism 

Let us consider a generic real signal described by the func- 
tion h(x), which depends on N parameters {p a }ae[ijf]- The 
measured signal is y(x) = h(x) + n(x), where n(x) is a realiza- 
tion of the noise specified by a probability distribution. Let 
us assume that the noise is Gaussian, is statistically station- 
ary with respect to x, has zero mean value, (n(x)) = 0, where 
() represents an ensemble average, and has known variance, 
a(x) 2 = (n(x) 2 ). The parameter estimation errors for p a can 
then be calculated using the Cramer-Rao bound [33] 



(<^ a <^>><r-V 



(49) 
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where equality is approached for high S/N signals. Here T a b 
is the Fisher matrix defined by 



Tab = 



d a h{x) dbh jx) 
a 2 (x) 



dx, 



(50) 



where d a is the partial derivative with respect to the parame- 
ter p a . Note that a(x) here is defined as the noise per unit x. 
In eq. (50), x denotes time f for time-domain samples, or / 
for frequency-domain samples. The purpose of this work is to 
study how an arbitrary end of the observation, at x max (or ff be- 
low, for time samples) affects the resultant correlation errors 
(Sp a Spb), for a fixed start-of-observation at x m j n (or t\ below, 
for time samples). 

An important quantity for the evolution of parameter esti- 
mation errors is the signal-to-noise ratio, S/N, defined by 



h 2 (x) 
a 2 (x) 



dx. 



(51) 



For LISA, the noise varies with signal frequency. In this 
case, the Fisher matrix can be evaluated in Fourier space [33, 
34], 



r(/D.i = » 4 



/' 



m B akf) dbhjf) 



df 



(52) 



where h(f) is the Fourier transform of hit), the GW signal 
(28), d a is the partial derivative with respect to parameter p a , 
bars denote complex conjugation, and S 2 (f) is the one-sided 
spectral noise density (§ II A). 



B. Approximate solution 



We are only interested in estimating uncertainties for the 
Psiow variables (§ II A), which are determined exclusively by 
h m (t). Recall from eq. (29) that \h c (t) \ = A so that, for the 
Fourier transform [40], we have \h c (t)\ 2 = 4\h c (f)\ 2 (df /dt). 
Using these relationships, let us define the instantaneous rela- 
tive noise amplitude per unit time <r(t) as 



CT - 2(f) = /cVo(0] 



S 2 n [fo(t)] 
3V5 



d/o(0 



df 



A 2 [f Q (t)] 

S 2 n lMt)] 



(56) 



4 S 2 n [f Q (M z ,t)Y 

The last equality follows from the Newtonian waveform and 
frequency evolution, eqs. (23) and (1). We point out that the 
mass dependence is captured entirely by a(t) and does not 
appear anywhere else in what follows. 

By combining eqs. (54), (55), and (56), we arrive at 



T(tf)al 



7" 

Jt f 



^[d a h m it) d h h m (t)] 
a 2 (t) 



df. 



(57) 



Equation (57) is the special case of (52), where the carrier 
signal-to-noise ratio and modulation, h m , are conveniently iso- 
lated. 

We are now ready to make use of the harmonic mode de- 
composition. Substituting eq. (31) into (57) gives 

T(tt) ab = ^\ J2 &8h8b8ji p hrh(tt)\, (58) 



7l ,72= 



where 



-df. 



(59) 



We seek an alternative equivalent form of eq. (52) specific 
to GW inspirals for which, as in eq. (28), the high frequency 
carrier signal is decoupled from the slow modulation. In case 
of SMBH inspirals, with a high frequency signal /z c (f) chang- 
ing its frequency slowly as /o(f) given in eq. (1), and fur- 
ther modulated by a slowly varying function h m (t) as given 
in eq. (28), the integral in eq. (52) can be evaluated in the sta- 
tionary phase approximation, by substituting 



h(f) = h m [t (f)]xh c (f), 



(53) 



where h c (f) is the Fourier transform of the carrier signal and 
h m [to(f)] is the modulating function evaluated at the time 
when the carrier frequency is /. This can be converted to the 
time domain, by simply changing the integration variable to 
t = to(f) using the frequency evolution in eq. (2): 



T( ti ) ab = 3iU 



f 

Jtf 



dghjt) d h h(t) 



d/o(0 



df 



df 



and 



h(t) = h m (f)x h c [f (t)]. 



(54) 



(55) 



The function Pj{k) is shown in Figure 1 for j = 0, to- 
gether with real and imaginary parts for the j = 4 case, for 
m\=nii = 10 6 M Q at z = 1. Since the accumulated signal-to- 
noise ratio is S/N = Po(t), the figure shows that the instanta- 
neous signal-to-noise ratio is [d/dt](S/N) = [d/dt]P Q (t) w f" 2 . 
The extrapolated signal-to-noise blows up at "merger" (t = 0). 
Data analysis for such a non-stationary signal-to-noise ratio 
evolution has several interesting implications, which we study 
further with toy models in Appendix A. We find that, for such 
a signal-to-noise ratio evolution, specific combinations of pa- 
rameters can always be measured to very high accuracy. 

The time dependence in eq. (59) couples only to the combi- 
nation j = ji—ji. This allows us to rearrange the double sum 
on (j'i , 72) and evaluate one of them independent of time: 

r(pWfU = 3? j X^Miow(0))UP;fe) j , (60) 
where 

8 

[^■(Psiow(0))Lfe = ^2 d a g j+f (p siov ,(0)) d b g j(p slov ,(0)). (61) 
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FIG. 1: The time evolution of the fundamental functions Pj{ti), used 
to construct the Fisher matrix to forecast localization errors by LISA. 
The dependence of the Fisher matrix on the look-back time tf is ob- 
tained from the 9 fundamental functions with < j < 8. The curves 
show Po(ff), as well as the real and imaginary parts of Pi,(tf), for 
mi = ni2 = 10 6 Mq and z = 1. Thin dotted lines represent negative 
values. Note that Poitf) = (S/N) , which is the simple scaling of 
inverse squared errors, neglecting correlations. The signal-to-noise 
ratio scales steeply, approximately as S/N oc tf 1 . The curve Pn(tt) il- 
lustrates how all the other similar Pj(tt) functions vary, with a relative 
number \ j\ of oscillations, and P-j(tt) = Pj(U). 



Our parameters in the correlation matrix are p a = 
(dL,0pf,4>N,0NL,4>NL) since we assume that the other param- 
eters, i.e. {A4,,?7z,</>isco,fmerger,a,7(0)}, are known from the 
high frequency carrier signal (§ II C). It is straightforward to 
compute the derivatives of gjid^, il) using eq. (42) for all pa- 
rameters p a , except <j>jf. The (j)^ dependence in gj in eq. (42) 
is implicit in a = a(E(0),$(0),(p N ) and 7(0) = 7($(0),</> JV ) 
(see § II A). Since we assume that S(0) and $(0) are mea- 
sured to very high precision with the high frequency carrier 
(§ II C), we can use the chain rule to express the <p N derivative 
as d l p N gj = d a gj-d^g j . 

Up to this point we did not make use of the fact that the 
LISA signal is equivalent to two orthogonal arm interferom- 
eters rotated by Aa = ir/4 with respect to each other. To ac- 
count for both data-streams being measured simultaneously, 
the Fisher matrix is written as the sum of the two Fisher matri- 
ces corresponding to each individual interferometer. Writing 
out only the a dependence, we have T^(a) = Y a b(a)+Y a b(a- 
7r/4). Finally, according to eq. (49), the parameter error co- 
variance matrix is the inverse of this total Fisher matrix: 

(SpaSpb) >[Tab(dh,0 N ,6 N L,<j>NL,a,j(0)) 

+ Y ah (d L , 9 N ,9 NL , (f> NL ,a-w/4, 7(0))]"' . (62) 

Equation (62) along with (60) is our final expression, describ- 
ing the time evolution of parameter estimation uncertainties. 
We note that after combining both data-streams, the matrices 
[.^/(PsiowCO))]^ for 4 < \ j\ < 8 modes vanish exactly for all 



Pslow(O). 

Let us emphasize the most important features of eq. (60): 

• The parameter dependence is separated from the time 
dependence. The Fisher matrix, Y a b, is written as a 
linear combination of matrices !Fj(p s i ov ,) weighted by 
the scalars Pj(t), where ^}(p s iow) is independent of time 
and Pj(t) is independent of the parameters p a . Evaluat- 
ing ^(psiow) requires the computation of the parameter 
derivatives d a gj- 

• The evaluation of all integrals Pj(tf n ) for different n = 
l,2,...,N tf can be done in the same amount of time as 
needed for a single integration, since the tf dependence 
enters only in the integration bound in eq. (59), 

• Large Monte Carlo (MC) simulations can easily be per- 
formed since the time evolution is given by a small 
number of functions, Pj(t), which can be calculated a 
priori and pre-saved. No integrations at all are neces- 
sary during the MC simulation for calculating distribu- 
tions of correlation matrices. 

In § V below, we estimate the improvement in the computa- 
tion time provided by the HMD method for calculating distri- 
butions of parameter errors and their time evolution. 



C. Generalization to the exact PN signal 

So far, we only considered the simplest case, assuming 
no cross-correlation between parameters p s i ow and (pf as t,.PspinX 
for a restricted post-Newtonian waveform. Moreover, we as- 
sumed the Doppler-phase (24) to be negligible. Including 
cross-correlations and the Doppler phase would allow us to 
examine the range of validity of our approximations, and it 
would allow us to extend computations to the lower compo- 
nent mass regime, m < 10 5 M Q , where the Doppler phase be- 
comes important. Furthermore, including spin precession ef- 
fects can modify angular localization errors by a factor of ~ 3, 
at least for the final errors at tf w / ISC0 [17, 22]. While we con- 
tinue to use our initial set of approximations in later sections, 
we outline here how the HMD formalism could be used to 
decouple the time-dependence from the angular parameter- 
dependence, even in the case of the most general (arbitrary 
order) restricted post-Newtonian waveform. Source sky po- 
sition angles (Q Nl cf) N ), detector angles at ISCO (a, 7(0)) and 
luminosity distance (<f L ) can all be decoupled even if spin- 
orbit and spin-spin precessions are included in the waveform, 
which is potentially a great advantage over the traditional cal- 
culation methods (see § V for a detailed discussion). 

Consider a general restricted post-Newtonian signal given 
by eq. (22), for which we substitute the harmonic mode ex- 
pansion [41], 



h hll (f) = _4^^ Ln /- 7 / 6 £ .'l^e'(/)+VD+0Gw] j 

j=-4 



(63) 
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where A = A(pf dS i) is the amplitude in the frequency-domain 
(eq. 68 in ref. [17]), gj = gj(p s i ov/ ) is the modulation ampli- 
tude in eq. (42), p D = </> D (Psiow,/) = c D (p dovj )f is the Doppler 
phase (see eq. 26 above), c D (p s i 0W ) = 27r/~' sin0 w cos7), 
(pew = 4>Gw(p',f) is the GW phase (20) (e.g. eq. 3.2 in 
ref. [25]), 



Equation (70) can now be substituted in the Fisher matrix 
in eq. (52). We get 



where 



r(ffW = r<™te)+i^te)+r$fe) 



CD/ 



<2), 



(71) 



*Gwfe/) = ^^ W 
«=0 



(Pfast,Pspin)«^ W (/), 



(64) 



and time-frequency relationships t (f) = t(p; /) can be written 
as (e.g. eq. 3.3 in ref. [25]) 



t(p;f) = ^cl{pf. isU p, via )ul(f). 



(65) 



n=0 



Here, the c GW and coefficients are frequency- 
independent, while m gw (/) and u]j(f) are parameter- 
independent. They correspond to the various post-Newtonian 
terms in the post-Newtonian expansion, and N corresponds to 
the highest order term. The frequency functions are very sim- 
ple powers of /, i.e. m gw (/) = /H5/3) and u T n {f) = /"-( 8 / 3 >. 
(Note that neither the c GW and c\ coefficients nor the m gw (/) 
and u^(f) functions are complex. Every term in eq. (63) is 
real except for the gj = gj(p s i 0VJ ) coefficients and the complex 
argument.) 

Equations (63-65) can be combined into 

/</) UI = £Af(p)r 7/ v*^/\ (66) 

j=-4 



where 



and 



Ay\p) = ^(pf as t)^' n (p s iow), (67) 



*,(?;/) = c D /+^E c ^»(/) + E c » Wii » GW W( 68 ) 

n=0 n=0 

2N+1 

= ^cjntttf), (69) 

k=0 

where in the last step we introduced u = {/, w T , w GW } and Cj = 
{c D , jw©c T ,c GW } to collect all /-functions and coefficients in 
one vector and one matrix, respectively. 

To compute the Fisher matrix, we need to obtain the partial 
derivatives of the signal with respect to the parameters: 

4 / 2N+1 \ 

~ h Af) = E A J> a + ' E -'•«.-»<</) r 7/ v*; (/) , (70) 

j=-4 \ k=0 ) 

where commas in indices denote partial derivatives with re- 
spect to the parameter following the index. Note, that the pa- 
rameter index a spans all parameters p s i ow > Pfast> an d ftpin, and 
the Fisher matrix accounts for correlations between these pa- 
rameters. 



(72) 



. h J2=-4 

4 2N+1 



7l,72=-4 k=0 



2N+1 



-i E A h A h,bY C h^ P hhk^ f > 



(73) 



h j'2=-4 k=0 

4 2N+1 



Tf b (tf)=3tl- E A h A h E C hkuaC hkl , b Pf hklk2 {k)\, 
{ jij2=-4 ki,k 2 =0 J 

(74) 

and where 



ffCO f-7/3 J{ji-ji)uj 9 t(f) 

CtW= 4r r7/3 ""^" d/. <75, 



/mi„ 5 «(/) 

/w /- 7/ x (/)«*2 (f)^-^"^ 

fmn 



s 2 „(f) 



df, 



are the frequency dependent terms. 

Equation (71) is our final result, where the localization 
parameters (i.e. angles and distance) are decoupled from 
all other parameters (i.e. masses, spins, reference time and 
phase at ISCO). The equation explicitly shows that, contrary 
to the traditional methods usually adopted for Monte Carlo 
computations of random binary orientations and sky positions 
[5, 7, 8, 13, 17, 19-22], the localization of a LISA inspiral 
event and its time-dependence can be explored without the 
need to evaluate integrals for each realization of the fiducial 
angles. Note that the only approximation made to obtain eq. 
(71) was to neglect of spin-orbit and spin-spin precession in 
the general restricted post-Newtonian solution for the Fisher 
matrix. The time-dependence is given by the P(ff) functions 
in eq. (75) and the extrinsic parameter-dependence is given by 
the coefficients, A. The P(tf) functions in eq. (75) can be com- 
puted a priori, independently of the fiducial angles. Note that 
P(tf) depends implicitly on the parameters (pf as t,Ps) through 
t(f) in eq. (65), and its inverse /(f). Generally, there are at 
most (2/ max +l)x [l+(2/V+l) + (l/2)(2/V+l)(2/V+2)] such 
independent functions. 

From the general case, we can now deduce the special so- 
lution in eqs. (60) and (62) valid for a Newtonian evolution, 
no Doppler phase, and no cross-correlations between p s iow 
and (pfast,Pspin)- This approximation simply corresponds to 
the first term r^(ff) in eq. (71), where in eq. (75) the t(f) 
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function is computed using the Newtonian formula fo(/) given 
by eq. (2). Note also, that the next term in eq. (71), T^fe), 
corresponds to the cross-correlation of the amplitude modu- 
lation with the "high frequency carrier signal" (i.e. Doppler 

(2) 

phase and GW phase). The last term, T a ^(tf), corresponds to 
the cross-correlations among parameters in the high frequency 
carrier. 

Finally, we briefly consider extensions to include spin-orbit 
and spin-spin precessions in the signal. Let us refer to the an- 
gular momentum angles as ph(t) = (Stuff), <j>m(t)), which are 
now time-dependent. As we briefly show next, in the case of 
spin precession, the P 01 ' 2 (/f) time-dependent integrals loose 
the convenient property of being independent of p^, but never- 
theless, the parameters describing the sky position and detec- 
tor orientation are still time-frequency independent and they 
are decoupled in this prescription. Indeed, an extra precession 
phase e\p[iS P (p L , p spin ,t)] has to be included in eq. (63) and 
gf (Psiow). We now have ^ vec ( P ;f) = * j(p;f)+8 ? (p uP ^ m ;t) 
instead of eq. (68). Thus, when taking the derivatives of 
the signal in eq. (70), we will get additional terms propor- 
tional to A I : 11 (p)d a 5p ta and the original A J fl terms will have 

time variation due to the p^ dependence of g*' n . Finally, after 
these modifications, the Fisher matrix will be similar to that 
in eq. (71), except now the p^ terms cannot be moved out- 
side of the frequency-integral but have to be attached to the 
time-varying P°' 1,2 (ff) part [42]. The main advantage we re- 
tain is therefore that the source position and detector angles 
(9N,4>N,a,-f(0)) will still only be included in the coefficients 
A 1 j u (p) and c^ a and the time-evolution can be still be com- 
puted independently of these parameters. In § V C we show 
that this indeed reduces computation times by a large factor 
relative to the traditional methods (e.g. [17, 22]). 

We leave numerical implementations and explorations of 
parameter distributions and their time-dependence, in this 
case of a general inspiral waveform, to future work. 



V. COMPUTATION TIME 

One of the great advantages of introducing the HMD 
method is the reduction in the computational time needed to 
evaluate error distributions for the parameters which deter- 
mine how efficiently LISA can localize SMBH binary inspiral 
events: sky position, angular momentum orientation and fi- 
nal detector orientation. In general, this is a computationally 
very demanding task because of the large dimensionality of 
the angular parameter space. Mapping the structure of the 
distribution of correlations in the parameter space of mock 
LISA measurements requires vast Monte Carlo simulations, 
which are presently limited by computational resources. Cur- 
rently, only a small portion of this space has been explored 
[7, 8, 13, 17, 20-22]. In this context, it is desirable to tune 
methods to the specific problem at hand. The HMD method 
described above is specifically constructed to exploit the struc- 
ture of LISA inspiral signals. 



A. Approximate solution 

The computational time for parameter space exploration, 
using the HMD method with the approximations described in 
§ IV B, is significantly reduced for the following reasons. The 
standard approach for estimating parameter errors requires 
the evaluation of an N p x N p symmetric Fisher matrix, where 
each matrix element is an integral over the range spanned by 
the GW frequency during inspiral (see Refs. [33, 34]; and 
§ IV). Here, N p is the number of parameters describing the 
signal. The number of independent elements in a symmetric 
matrix is (1 /2)N P (N P + 1). Let us assume that the evaluation 
of a single integral requires to compute the waveform at N nt 
separate instances. The computation of one integral is suffi- 
cient also to trace the time evolution at N nt different ff val- 
ues, if one uses a single trapezoidal integral in frequency from 
/isco to / m i n and stores results at each intermediate value of 
/. Since the time evolution of the frequency is known inde- 
pendently of the angles, we can already get the integral for 
N nt different ff values. For randomly chosen fiducial angles 
in a MC simulation of size Nmc requires the calculation to 
be repeated Nmc times. To evaluate the evolution of param- 
eter errors as a function of ff for N k different instances re- 
quires N k computations. Therefore, the standard method costs 
Ar standard = (\ /2)N p (N p +\)N mt N M c computational time units. 

In contrast, with our proposed HMD method (§ IV B), 
the p s iow parameters are decoupled from the pf. ist parameters 
(§ II C) and from time. The N p x N p Fisher matrix can be 
split into two smaller matrices, with N p \ x N p i and N p q x Npa 
components, where N p = N p \ +N p q. The N p \ x N p \ matrix de- 
termines the angular errors (which are deduced from the am- 
plitude modulation of the signal), while the other matrix deter- 
mines the remaining parameters (e.g. masses, phase and time 
at ISCO, using the high frequency carrier only). Since we are 
only interested here in parameters relating to the localization 
of the source, p s iow, it is sufficient for us to consider the corre- 
sponding N p \ x N p i matrix only. Since it is symmetric, it has 
only (1 /2)Np\{N p \ + 1) independent elements, but it turns out 
that the computation of only N p i elements is sufficient for a 
single harmonic (we need only the N p \ derivatives of the gj 
functions, see eq. [58]). Using 2Nj+ 1 harmonic modes and a 
MC simulation with Nmc random choices of fiducial parame- 
ters costs N p \{2Nj+ \)Nmc time units. In this method, the time 
dependence is decoupled, so that parameter dependencies can 
be taken outside of the integral (see eq. [58] and Tj(p & \ ovl ) in 
eq. [60]). The MC sampling can therefore be evaluated inde- 
pendently of time. The time evolution of the signal for each 
harmonic is known independently of the angles, by construc- 
tion (Py(ff), eq. [59]), and this integration for each component 
can be evaluated a priori, independently of the fiducial param- 
eter values. For each such mode, we would like to evaluate a 
number N k of integrals. Fortunately, since the different inte- 
grals differ only via the lower integration bound in the time 
domain, all integrals can be obtained during a computation of 
the integral with the largest time domain, ff = fisco- There- 
fore for a total of {2Nj+ 1) modes, building the time-evolution 
functions Pj(ti) takes of order (2Nj + l)N tl time units. This is 
generally much faster than building the time-independent co- 



13 



efficients ^/(Psiow)' In summary, with the HMD method, one 
only needs Ar HMD = N pX (2Nj + l)N MC + (2Nj + l)N tf units of 
time. 

Comparing methods, we find that the computational re- 
quirements of the HMD method is lower by a factor of 



A7 HMD m _ N pl (2Nj+l)N M c + (2Nj+ l)N t! 



a rj.no spin 
standard 



(l/2)N p (N p +Win t N MC 



(76) 



Recall from § II A that the number of parameters for a no- 
spin case is N p q = 4, N p i = 5, so that N p = N p q+N p \ = 9. 
Choosing N int = 10 3 , N U c = 10 4 , N t{ = 100, and Nj = 4 for 
the other parameters in eq. (76) as a representative example, 
the gain in computational efficiency is Ar st a n dard/A7HMD = 
(5 x 10 8 )/(5 x 10 5 ) w N int = 1000. Moreover, the Fisher ma- 
trix is much smaller, 5x5, which offers an important further 
advantage when performing the inversion to obtain the error 
covariance matrix. Using the HMD method, the inversion of 
the Fisher matrix is computationally less expensive than gen- 
erating the matrix. Note that the second term in Ar^p" 1 , 
corresponding to the Pj(ti) functions, is negligible in this case 
and the computation time is dominated by constructing the 
coefficient matrices. A calculation of the representative MC 
example above, with a non-optimized implementation of the 
HMD method, takes less than a minute on a regular worksta- 
tion. 

The case for substantial improvement with the HMD 
method becomes even more compelling when additional pa- 
rameters (spins and higher order PN terms) are included, as 
we discuss next. 



B. Post-Newtonian Signal without Spin Precession 

We now consider the general HMD method outlined in 
§ IV C, with N S pi n = 6 spin components. The spin parame- 
ters can be grouped as Nsm = 2 independent spin magnitudes 
and Nsa = 4 independent spin angles. Since the spins can 
be oriented arbitrarily, the spin angular parameters have to 
be randomly chosen, in addition to the other angular param- 
eters in any Monte Carlo computation. This enlargement of 
the parameter space of random parameters greatly increases 
the computational cost, both for the standard method and the 
HMD method. However, we show next that the incremen- 
tal cost is much less severe for the HMD method. The HMD 
method should be considered in future work aimed at comput- 
ing time-dependent parameter estimation errors in the general 
case with spins. Here, we neglect the effects of spin preces- 
sion, so that the angles ft = (6N,fa,0NL,faL,ce,j(O)) are de- 
coupled and, since the signal does not depend on fa, there are 
only Nq = 5 independent (spin-unrelated) angles. 

The larger the parameter space, the larger the sample size 
must be in a Monte Carlo computation. Let us assume that the 



AT 

1 V . 



(i) 



where N^ c is the 



sample size is chosen to be Nmc - y>MC 
effective number of samples for a single parameter, and d = 
Na when spanning the _!-space only, d = Nsa when spanning 
the spin-angle space only, and d = Nn+ Nsa when spanning 
both. 



To compute the time-independent matrices, we have to 
evaluate N p (2Nj + 1) independent A j A coefficients for the full 
f2-space and we have to compute the N p (2Nj+ l)(2iV fw + 1) 
independent cjk .„ matrices over a d = No + Nsa dimensional 
parameter space. (Here No = 2 denotes the number of param- 
eters on which the Doppler phase depends, (8 N ,-f) in eq. [26], 
using the fact that both cgw and c T also depend on all spin 
angles in eqs. [64,65].) In § IV C we have shown that there 
are (2Nj+ 1) x (2N pn +7Npn + 6) independent integrals, where 
Npm is the number of terms in the post-Newtonian expansion 
plus the Doppler phase. We have to compute these integrals 
for all spin angle orientations. Therefore, the computational 
cost scales as 



a 7.110 spin prec _ » r 
HMD ~ ly P 



0) 



(2N J+ l)(i. MC j 
+(2Nj+l)(2N 2 PN +lNp N + 6)N t[ 

Nd+Nsa 



Nsa 



+N p (2Nj+1)(2N pn +1)(n£> c ) . (77) 

For the standard method, reiterating the argument in § VA, 
we get 



a 7.110 spin prec 
standard 



1 / m \Na+Nsi 

= ^N p (N p +l)N int [N^ 1 



"mcJ 



(78) 



where now N p = N pl +N p0 +7V S p in = 15. Taking N^ Cl = 10, 
N PN = 4, Nj = 4, N D = 2, N n = 5, N SA = 4, N mt = 10 3 , 
N t[ = 100 as a representative example, we find that the 
HMD method is computationally less expensive by a factor 

A Cnd P ard Pre 7 A7 HMD in *"* = d X 10 14 )/(2 X 10 9 ) = 7 X 10 4 , 

as compared to the standard method. 

Note that the first term in eq. (77) corresponds to our orig- 
inal approximations in § IV B, i.e. no cross-correlations be- 
tween p S pi n and Psiow This approximation indeed leads to 
much faster computations, since Ar s "° n ^ prec /Ar^ D DXlm '' te = 
9 x 10 6 for the same representative example. 



C. Post-Newtonian Signal with Spin Precession 

Accounting for spin precession, the Nl = 2 angular momen- 
tum angles, (6 fa£), and the Nsa s P m angles are now chang- 
ing with time. In this case, one has to solve a differential 
equation for the evolution of these angles for each individual 
Monte Carlo realization. We assume that this can be com- 
puted independently of the Fisher matrices and that it would 
take NoEN t[ computation time units to evaluate, at each of the 
N tf time instances, for each initial set of angles. 

The HMD method costs 



Ar™° =N p (2N J+ l)(N^ c 



n N d +Nsa 



Nl+Nsa 



+N p (2Nj+l)(2N PN +2)(N^ c j 
+(2Nj+ l)(2N 2 PN +7N PN + 6)N, (_#>,)' 

/ n s\ Nl+Nsa 

+N DE N t! [N^cj , (79) 
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where the first term involves constructing the Aj M coefficient 
matrices, the second term involves constructing the c^ a coef- 
ficient matrices, the third term is for computing all three time- 
evolution quantities P(tf) in eq. (75), and the fourth term is for 
solving the precession equations. 

In the standard method, we need first to solve the precession 
evolution differential equation and then construct the Fisher 
matrix [22]. Following the assumptions made above, we esti- 
mate a cost 

Ar™ = N DE N tl (a£») ' 

+-N p (N p +l)N mt [N^ c j N tr (80) 

Using Nl = 2 and all other parameters as in § V B, we find 
that the HMD method is computationally more efficient by a 
factor Ar"d7Ar™ c = (2 x 10 14 )/(6 x 10 10 ) = 2 x 10 3 , 
as compared to the standard method. Since the (9nl,4'nl) 
subspace could no longer be decoupled, the efficiency of the 
HMD method relative to the traditional method lost a factor 
of 30, as compared to the no spin precession case in § VB. 
Nevertheless, the computational advantage remains very sub- 
stantial. 



VI. RESULTS 

Having described the HMD formalism in detail, we now 
apply it to build MC simulations aimed at studying how RMS 
source localization errors [43] evolve as a function of look- 
back time, ff, before merger. The low computational cost of 
the HMD method allows us to survey simultaneously the de- 
pendencies on source sky position, SMBH masses and red- 
shifts. We carry out MC calculations with 3 x 10 3 random 
samples for the angles co&9n,cos9nl,4'nl,C(,7(0)- Several 
thousands values of M and z are considered, in the range 
10 5 < M/M < 10 8 and 0.1 < z < 7. In addition, we ran 
specific MC calculations to study possible systematic effects 
with respect to the source sky position, by fixing 9 N and (f> N 
(on a grid of several hundred values) and varying all the other 
relevant angles. 

In all of our computations, we calculate the er- 
ror covariance matrix for c/l, #a<, 9nl,<I>nl- Following 
Lang & Hughes [22], we calculated the major and minor axes 
of the 2D sky position un certa inty ellipsoid, 2a and 2b, and 
the equivalent diameter, \J\ab. 

We have verified our HMD implementation and the gen- 
eral validity of our assumptions by comparing our results at 
ISCO with those of Lang & Hughes [22] (for m\ = mi = 
10 5 , 10 6 , 1O 7 M at z = 1 and m x = m 2 = 10 5 , 1O 6 M at z = 3, 
in the no spin precession case). Depending on SMBH masses 
and redshifts, we found agreement at the 5-30% level for the 
mean errors on the luminosity distance, major axis, and minor 
axis. The small discrepancies may be due to differences in the 
set of assumptions made. Lang & Hughes [22] account for 
the small cross-correlations between the p slow and {pf ast ,Pspin} 
parameters and they choose f; to be uniformly distributed be- 
tween merger time and LISA'S mission lifetime. Recently, 



Lang & Hughes reported angular errors that are a factor of 
2-3 lower [35], which are inconsistent with our results at this 
level. Nevertheless, these discrepancies are still small relative 
to the typical width of error distributions or to the systematic 
variations of mean errors with tj, M, and z (from a factor of 
few to orders of magnitudes, see Fig. 2 below). This success- 
ful comparison justifies the use of the HMD method to study 
the dependence of localization errors on look-back time. 



A. Time dependence of source localization errors 

We calculate the variation with look-back time, ff, of the 
distribution of marginalized parameter errors for a range of 
values of M,z, 9 N ,9 NL ,(j) NL ,a,j(0). Figure 2 shows results for 
random angles and rti\=m 2 = 10 6 M Q , at z = 1. 

The top panel shows the luminosity distance error, dd^, 
while the bottom panel describes the equivalent diameter, 
2^/ab, of the sky position error ellipsoid with minor and ma- 
jor axes a and b. The figure displays results for three separate 
cumulative probability distribution levels, 90%, 50%, 10%, so 
that 10% refers to the best 10% of all events, as sampled by 
the random distribution of angular parameters. The evolution 
of errors scales steeply with look-back time for ff > 40 days. 
In this regime, the improvement of errors is proportional to 
(S/N)~ l . For smaller look-back times, errors stop improving 
in the "worst" (90% level) case, improve with a much shal- 
lower slope than (S/NT 1 for the "typical" (50% level) case, 
and keep improving close to the (S/N)' 1 scaling in the "best" 
case (10% level among the realizations of fiducial angular pa- 
rameters). Although Figure 2 shows only the equivalent di- 
ameter of the 2D sky localization error ellipsoid, we have also 
computed the evolution of the distribution of the minor and 
major axes. We find that amb w ^/ab initially (i.e. the ellip- 
soid is circular), but the geometry changes significantly during 
the last two weeks to merger. For example, in the typical case, 
the major axis a stops improving at late times, while the minor 
axis a maintains a steep evolution. Therefore the eccentricity 
of the 2D angular error ellipsoid changes quickly with look- 
back time. This is important because large eccentricities can 
play a role in assessing observational strategies for EM coun- 
terpart searches [18]. 

To map possible systematic effects with respect to source 
sky position, we carried out MC computations with ran- 
dom (cos9 NL ,4> NL ,a) angles (sample size Nmc = 3 x 10 3 ) but 
fixed source sky latitude and longitude relative to the detec- 
tor (9 N ,-f), for mi =/W2 = 1O 6 M and z = 1. We find no sys- 
tematic trends with sky position for Sd^, for any value of the 
look-back time, ff . Neither do we find systematic trends with 
sky position for the distributions of minor and major axes of 
the angular ellipsoid, for any value of the look-back time, ff, 
as long as 9 N is not along the equator. The case of equatorial 
sources, with 9 N w 90° and a short look-back time ff before 
merger, is the only nontrivial one we have identified. In that 
case, we find a minor systematic trend with 7 longitude. The 
error distributions shift periodically up and down, relative to 
the average, when changing 7 from to 2ir. 

In addition, to map dependencies with mass-redshift-look- 
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FIG. 2: Evolution with pre-ISCO look-back time, tf, of LISA source 
localization errors, for M = 2 x 1O 6 M and z = 1. The top panel 
shows luminosity distance errors and the bottom panel shows sky 
position angular errors (equivalent diameter, 2^/ab, of the error el- 
lipsoid). Best, typical, and worst cases for random orientation events 
represent the 10%, 50%, and 90% levels of cumulative error distribu- 
tions, respectively. Errors for worst case events effectively stop im- 
proving at a finite time before ISCO, even though the signal-to-noise 
ratio accumulates quickly at late times. Errors for best case events 
(especially the minor axis) follow the signal-to-noise ratio until the 
final few hours before merger. 



back time of localization errors, we carried out MC compu- 
tations with arbitrary (cos6s,cos8 NL ,(f) NL ,a,-f) angles, with 
sample size Nmc = 3 x 10 3 , for several thousand pairs of 
(M,z) values. We find that the evolution with look-back 
time of error distributions depends sensitively, and in a 
complicated way, on the mass-redshift parameters. Gener- 
ally, localization errors increase with redshift. Firstly, the 
S/N is approximately proportional to the instantaneous value 
cT(h S co)^V 3/4 l(l+z)M] 5 / 4 /ddz)S n (f lsco r 1 (eq. [56]) and, 
secondly, the beginning-of-observation time scales as t\ oc 
rf l \{\ +z)M]" 5 / 3 (eq. [2]). For (1 + z)M < 4 x 1O 6 M , the to- 
tal observation time can exceed one year and the second effect 



is unimportant. We further describe mass-redshift dependen- 
cies below, in § VI B, in relation to advance warning times for 
targeted electromagnetic counterpart searches. 

The results on localization errors from our extensive explo- 
ration of the parameter space of potential LISA sources can 
be summarized as follows: 

1 . Probability distributions 

• The error distributions for 5d L , 2a, and 2b all have 
long tails: 1 %— 99% cumulative probability levels 
are separated by factors of ~ 100, while the 10%- 
90% levels are separated by factors of <~ 10. 

• The <5c/l distribution is skewed, with a median 
closer to the best case, a median smaller than 
the mean, even on a logarithmic scale. On the 
other hand, sky localization error distributions are 
roughly symmetric on a logarithmic scale. 

2. Fiducial parameter dependencies 

• The <5fl?L errors are roughly independent of fiducial 
angles throughout the observation. 

• For non-equatorial sources, the distribution of sky 
localization errors, (2a, 2b), is independent of sky 
position, i.e. the distribution does not have a sys- 
tematic dependence on 6s and 7 = $-4>n (for ran- 
dom a, 9 NL ,(j> NL ). 

• There is a small systematic trend with 7 for equa- 
torial sources. 

• There is a complicated dependence of sky lo- 
calization errors on M, z, and look-back time 
t f . For (1 +z)(ry/0.25) 3 / 5 M < 4 x 10 6 , and long 
observation times, errors scale with (S/N)' 1 w 
[(1 +z)M]^ 4 d L (z)-'S„(f lsC o(M,z))- 1 fa(td, where 
fa(tf) is the ff-scaling shown in Fig 2. For larger 
redshifted masses, the scaling has a complicated 
structure in the M, z, tf space that we did not ana- 
lyze in detail (but see eq. (A10) in the Appendix 
for scalings in terms of t\ and ff.) 

3. Time dependence 

• Luminosity distance and sky localization errors 
roughly scale with (S/N)~ l until 2 weeks before 
ISCO. 

• For the luminosity distance Sd^ and the major 
axis 2a, there is little improvement within the last 
week before ISCO for the typical to worst cases 
(i.e. 50%-90% levels of cumulative error distri- 
butions). 

• For the minor axis 2b, only the worst case events 
stop improving within the last week. The typical 
to best cases continuously improve until the last 
hour. 

• The eccentricity of the sky localization error el- 
lipsoid changes with time during the first and last 
two weeks of observation. The eccentricities are 
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smaller in between these two time intervals. For a 
detailed discussion of the eccentricity and its im- 
pact on counterpart searches, see Ref. [18]. 

• For the luminosity distance <Wl, the relative width 
of error distributions does not change during ob- 
servation and variations in the difference between 
the 90% and 10% levels of the cumulative distri- 
butions do not exceed 10%, except for the initial 
weeks, when the distribution is much more spread 
out. 

• For the sky localization errors, the width of error 
distributions increases during the final two weeks 
of observation, by a factor ~ 2 for the major axis 
and a factor <~ 4 for the minor axis. 



B. Advance warning times for EM searches 

From the astronomical point of view, being able to identify 
with confidence, prior to merger, a small enough region in the 
sky where any prompt electromagnetic (EM) counterpart to 
a LISA inspiral event would be located, is of great interest. 
With sufficient "warning time," it would then be possible to 
trigger efficient searches for EM counterparts as the merger 
proceeds and during the most energetic coalescence phase. In 
particular, an efficient strategy to catch such a prompt EM 
counterpart would be to continuously monitor with a wide- 
field instrument a single field-of-view (FOV), through coales- 
cence and beyond. Astronomical strategies for EM counter- 
part searches are the focus of a second paper in this series 
[18]. 

Given an angular scale, #fov, corresponding to the hypo- 
thetical FOV of a specific astronomical instrument, it is thus 
of considerable interest to determine the value of the look- 
back time ff at which the major axis, minor axis or equivalent 
diameter of the sky localization error ellipsoid provided by 
LISA just reach the relevant 6* F ov scale. This would allow one 
to trigger an efficient search for EM counterparts, in a well 
defined region of the sky that can be monitored. We will here- 
after refer to this time as the advance warning time. Note that 
it is important to differentiate the sizes of the major and mi- 
nor axes of the angular error ellipsoid in this context because 
the eccentricity can be large, and thus important in assessing 
optimal strategies for EM counterpart searches [18]. 

For definiteness, we evaluate advance warning times for an- 
gular diameters 6>fov =1° and 3.57° here but generalizations 
to other 6*fov values are obviously possible. The choice of the 
latter figure is motivated by the lOdeg 2 FOV proposed for the 
future Large Synoptic Survey Telescope, or LSST [36]. Fig- 
ure 3 shows advance warning times for a fixed source redshift 
at z = 1 and various values of the total SMBH mass, M. Fig- 
ure 4 shows similar results for various source redshifts, at a 
fixed value of M = 2 x 1O 6 M . 

In each case, we consider equal mass SMBH binaries and 
a maximum observation time of 1 yr (or lower if set by the 
GW noise frequency wall at 0.03 mHz). Each panel in Figs. 3 
and 4 shows the values of advance warning times at which 
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FIG. 3: Advance warning times (in days) for equal mass binary in- 
spirals at z = 1, as a function of total mass, M (in solar units). Best, 
typical, and worst cases refer to 10%, 50%, and 90% levels of cu- 
mulative error distributions for random orientation events, as before. 
The advance warning times shown correspond to the values of look- 
back times when the equivalent diameter, 2\fab, of the error ellip- 
soid first reaches 1° (top panel) or an LSST-equivalent field-of-view 
(3.57°, bottom panel). In the top panel, the worst case events are 
not shown because angular errors are too large even at ISCO. For the 
largest mass SMBHs, the maximum observation time (and thus t\) is 
below one year. 



the equivalent diameter 2\fab of the localization error ellip- 
soid drops below the reference #pov value. For each case, 
we show results for cumulative error distribution levels of 
10%, 50%, and 90%, labeled "best", "typical," and "worst" 
cases, as before. Figure 3 shows that LISA can localize 
on the sky events at z = 1 to within an LSST FOV at least 
one month ahead of merger, for 50% of events with masses 
2 x 1O 5 M < M < 3 x 1O 6 M , and at least 4 days ahead 
of merger for 90% of events in the same mass range. Fig- 
ure 4 shows that advance warning times decrease with red- 
shift, leaving at least 1 day ahead of merger for 50% of events 
with M = 2 x 1O 6 M , as long as z < 1 for 6> FO v = 1° and as 
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FIG. 4: Same as Fig. 3, for a fixed total mass M = 2 x 1O 6 M but 
various values of the source redshift, z. 



long as z < 3 for an LSST FOV. For events with this mass 
scale and the LSST FOV, there is a 10% chance that a 1 day 
advance warning is possible up to z ~ 5-6. 

Figures 3 and 4 display advance warning times for single 
one dimensional slices of the full {M,z) space of potential 
LISA events. With the HMD method, however, it is possi- 
ble to explore the entire parameter space of SMBH inspirals 
by repeating the calculation on a dense grid of (M,z) val- 
ues. We construct a uniform grid in the (logM,z) plane, with 
Az = 0.1 and AlogM = 0.1, and perform MC computations 
with 3 x 10 3 randomly oriented angles for each grid element. 
As a result, we obtain a complete description of the time evo- 
lution of sky localization errors in the large parameter space 
of potential LISA sources. Figure 5 displays advance warning 
time contours from this extensive MC calculation, for typical 
(50%) and best case (10%) events, adopting the LSST FOV as 
a reference. 

Advance warning time contours are logarithmically spaced, 
with solid-red contours every decade and a thick red line high- 
lighting the 10 day contour. Since advance warning times 
were computed on a finite mesh, contour levels for arbitrary M 




FIG. 5: Contours of advance warning times in the total mass (M) and 
redshift (z) plane with SMBH mass ratio m\/m,2 = 1. The contours 
trace the look-back times at which the equivalent radius {2\fab) of 
the localization error ellipsoid first reach an LSST-equivalent field- 
of-view (3.57°). The contours correspond to the 50% (top) and 
10% (bottom) level of cumulative distributions for random orienta- 
tion events. The contours are logarithmically spaced in days and 10 
days is highlighted with a thick line. 



and z values were obtained by interpolation. Our interpolated 
mesh is smooth if ff ~ 0.1 day, but it gets edgy for short ad- 
vance warning time approaching ISCO. Figure 5 shows that a 
10 day advance warning is possible with a unique LSST-type 
pointing for a large range of masses and source redshifts, up 
to M ~ 3 x 1O 7 M and z ~ 1.9. The bottom shows how far 
the advance warning concept can be stretched, by focusing on 
the 10% best cases of random orientation events. In this case 
a 10 day advance warning is possible up to z ~ 3 for masses 
around M ~ 10 6 M Q . Note that, in both cases, allowing for a 
warning of just one day would extend considerably the range 
of masses and redshifts for which a unique LSST-type point- 
ing is sufficient. 

These results can also be generalized to unequal-mass 
SMBH binaries. At fixed total mass, M, an unequal-mass bi- 
nary has an instantaneous signal-to-noise ratio that is reduced 
because of a lower rj value, but it also has a total observation 
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time that is potentially longer. Localization errors for unequal- 
mass inspiral events with total observation times longer than 
a month (i.e. with rfc ' 2 5 5 (l +z)M < 1 .8 x 10 7 M Q ) are degraded 
relative to the equal-mass cases discussed so far. For larger 
total mass, however, the worsening of errors is mitigated, or 
even reverted, relative to the equal mass case, thanks to the 
longer observation time. The error ellipsoid also becomes less 
eccentric thanks to this additional observation time. Figure 6 
summarizes results on advance warning times from the same 
MC computations as in Fig. 5, but this time for unequal-mass 
SMBH binaries with mass ratio m\/m.2 = 10. Despite a sys- 
tematic degradation in advance warning times (especially no- 
ticeable at low M values), the main effect of introducing a 
mass ratio m\/m2 = 10 is to shift advance warning time con- 
tours to somewhat larger values of total mass, M. Our main 
conclusions on advance warning times are not very strongly 
affected by the inequality of mass components in the popula- 
tion of SMBH binaries considered. 

Finally, it is important to understand how sensitive the re- 
sults are to the LISA detector characteristics. In particular, we 
examined how advance warning times are affected by increas- 
ing the minimum frequency noise wall or by loosing one of 
the arms of the 3-arm constellation. Figure 7 displays results 



for / m in =10 4 Hz, for m\/m.2= 1 and nti/mi =10. Increas- 
ing /nun mostly reduces the total observation time for high 

—Q/"X _ < j'X 

mass inspirals (?; ~ f mi ^ M z ; see eq. [2]) and reduces the 
signal-to-noise ratio by a small factor. As a result, the ad- 
vance warning time contours primarily shift in the (M,z) plane 
in the direction of smaller total masses by a factor of ~ 7, 
and secondly shift moderately (30-50%) to smaller redshifts. 
Loosing one LISA arm (i.e. using only one of the two interfer- 
ometers) most importantly removes the ability of the second 
datastream to break correlations in localization errors and also 
reduces the signal-to-noise by a small factor. As a result errors 
do not improve much during the last ~ 10 days before merger. 
Compared to the case with two interferometers, contours rep- 
resenting an advance warning of less than 10 days are shifted 
to significantly smaller z (especially for the minor axis of the 
sky localization ellipsoid), close to the 10-day contour, but 
warning times beyond ~ 10 days worsen only moderately. We 
conclude that even if / m j n = 0.1 mHz or if only one of the two 
interferometers is used, LISA still admits 10-day advance lo- 
calizations for a broad range of masses and redshifts, between 
10 5 <M<2x 10 6 andz< 1. 
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VII. DISCUSSION 

We have introduced a novel technique, the HMD method, to 
compute time-dependent GW inspiral signals for LISA. The 
method relies on the fact that LISA'S orbital motion induces 
a modulation on timescales that are long relative to the in- 
spiral GW frequency. Since this modulation is periodic, with 
a fundamental frequency of /©, it can be expanded in a dis- 
crete Fourier sum. In the HMD formalism, dependencies on 
sky position, orbital angular momentum orientation, and de- 
tector orientation in the LISA signal are inscribed in time- 
independent coefficients, while time-dependent basis func- 
tions are independent of these angles. This decomposition 
helps to reduce the computational cost of Monte Carlo sim- 
ulations exploring the time-dependence of source localization 
errors by orders of magnitude. 

Moreover, the HMD method can be used in conjunction 
with plausible approximations to further decrease the com- 
putational cost of explorations of the parameter space of lo- 
calization errors for LISA inspiral events. In our analysis, we 
identified two different characteristic frequency constituents 
of the signal: the high frequency restricted post-Newtonian 
GW inspiral waveform and the low frequency amplitude mod- 
ulation resulting from the detector's orbital motion. In the 
HMD method, these two components separate and parameters 
that depend only on the low frequency modulation (such as the 
source position and the orbital angular momentum angle) can 
be estimated independently of the other source parameters de- 
termined by the high frequency carrier signal. Our working 
assumption was that cross-correlations among these two sets 
of parameters must be much smaller than parameter correla- 
tions within either set. This hypothesis is valid very generally 
in the no spin limit for SMBHs, as shown by full Fisher matrix 
calculations without such approximations for general relativ- 
ity [7] and alternative theories of gravity [13]. 

In order to further examine the validity of our assumptions 
and the ultimate boundaries of our models, and to understand 
our results, we have constructed illustrative toy models that 
we now describe in some detail. These toy models show 
that the separation of parameters into various subsets associ- 
ated with different characteristic frequencies of the signal is a 
rather general property, which turns out to be an efficient way 
of reducing the computational cost of error estimations for the 
LISA problem. 

A. Simple toy models 

In this section, we discuss very simple toy models which 
capture the essence of the problem posed by the time- 
evolution of parameter error estimations. We then use these 
models to answer general questions on the LISA-specific pa- 
rameter estimation problem. 

Our harmonic decomposition technique is based on the sim- 
ple intuition that the angular information can be deduced from 
the slow periodic modulation of the high frequency GW wave- 
form. In § III, we have shown that modulation harmonics with 
frequencies larger than 4 / e vanish exactly. Here, we discuss 



the general properties of such a modulation. In the case of 
LISA, the high frequency carrier signal has an effective, cycle- 
averaged signal-to-noise ratio which monotonically increases 
with time as SMBH binaries approach merger. To mimic such 
events, we also assume in all of our toy models that the instan- 
taneous signal-to-noise ratio continuously improves through- 
out the observation. 

We seek answers to the following questions: 

1 . How do mean errors evolve during the final days of ob- 
servation? 

On the one hand, in standard angle-averaged treatments 
(e.g. [5, 13, 15]), an evolution of errors with the inverse 
of the signal-to-noise ratio is generally assumed. This 
would suggest a large improvement during the last day 
of inspiral. On the other hand, the slow modulation pic- 
ture suggests just the contrary: not much improvement 
is expected at late times when there is effectively very 
little modulation (Finn & Larson 2005, private commu- 
nication). 

2. Does the introduction of additional high frequency 
components in the signal have any effect on the esti- 
mations of low frequency parameters? 

In the GW context, it is of general interest to deter- 
mine under what circumstances additional high fre- 
quency signal components, such as higher order post- 
Newtonian corrections or spin-induced effects, remain 
decoupled from the determination of angular and dis- 
tance information based on the signal amplitude modu- 
lation. 

3. Are there combinations of signal parameters for which 
errors improve rapidly in the last days of observation? 
If so, what are these combinations? What determines 
how many such rapidly-improving combinations there 
will be? 

If the distance <f L correlates with the angles, then in 
principle the volume of the 3D error box can be much 
smaller than the product of the marginalized errors 
<5S1 x <5c/l would imply. Unfortunately, in practice, this 
is unlikely to help to reduce the number of false coun- 
terparts, because the Sz error will be dominated by weak 
lensing [15]. 

4. How does the width of parameter error distributions 
evolve with time? Are the best and worst cases ap- 
proaching the typical case prior to the final days of ob- 
servation? How do we expect the eccentricity of local- 
ization error ellipsoid to evolve with time for LISA? 

Here, we restrict our discussion to a brief summary of our 
findings and direct the reader to Appendix A for further details 
on these toy models. 

The parameter estimation uncertainties are defined by the 
correlation error matrix. For N p parameters, this defines an 
N p -dimensional error-ellipsoid in the jVp-dimensional param- 
eter space, where parameters are constrained at a given con- 
fidence level. Marginalized errors for a given parameter are 
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then related to the projection of this ellipsoid on the basis 
vector corresponding to that parameter. Since the principal 
axes of this error ellipsoid are generally not aligned with the 
original parameters, the marginalized errors can be substan- 
tial even if the volume of the error ellipsoid is close to zero. 
This happens if the ellipsoid is very "thin" but has a large size 
in at least one direction. Diagonal elements of the correlation 
matrix provide marginalized squared errors on the parameters, 
while eigenvalues provide squared errors along the principal 
axes. 

We consider three versions of toy signals to understand 
how a particular harmonic mode contributes to the time- 
dependence of parameter uncertainties and to find answers to 
Questions 1-4 above. We start with the simplest toy model 
and refine this model by adding more details and complexity 
in the successive models. In each case, we discuss general 
implications for the model under consideration. 



1. Basic toy model 

In our basic toy model, we assume that the true signal is 
comprised of a constant carrier signal, which is modulated by 
a single known-frequency cosine, /© : 

/z(f) = c + c 1 cos(27r/©f), (81) 

where cq and c\ are unknown parameters to be estimated. We 
assume that the noise level is rapidly decreasing during the ob- 
servation, mimicking the gradual increase in the instantaneous 
signal-to-noise ratio for LISA inspiral signals. The contradic- 
tory statements made in relation to Question 1 above can be 
explored with this model. We find that marginalized param- 
eter errors scale with the signal-to-noise ratio far away from 
merger (i.e. tf > O.l/^ 1 ) but they quickly converge to their fi- 
nal values at late times, even though the signal-to-noise ratio 
keeps accumulating. It is possible to derive analytical formu- 
lae for the evolution of parameter errors to fully characterize 
this behavior (see Appendix A). We find that, even though the 
error ellipse rapidly decreases in volume, as the inverse of the 
signal-to-noise ratio near merger, the error ellipse only shrinks 
along one of its dimensions, the semi-minor axis, so that a 
non-negligible residual uncertainty remains in the orthogonal 
subspace (e.g. along the semi-major axis). This residual un- 
certainty carries over to final marginalized errors for both pa- 
rameters. Therefore, this first toy model verifies the second 
option in relation to Question 1.: there is no late improvement 
because there is very little effective signal modulation, mak- 
ing the signal-to-noise argument largely irrelevant. However, 
we find below that this model does not carry some essential 
features of the LISA signal which modify somewhat our final 
answer to Question 1 (see final toy model below). 



2. Second toy model 

In our toy second model, we modify the single frequency 
signal by postulating two pairs of unknown amplitudes and 



phases for two different a priori known frequencies, satisfying 
fi ^> fx, which modulate an otherwise constant signal: 

h(t) = co + si sin(27r/if) + ci cos(2nf\t) 

+ Si sin(27r/2?) + CioCOs(27r/ 2 0. (82) 

The number of unknowns in this model is five: cq, s\ , c\ , Sio 
and cio are the coefficients of the functions 1, sin(27r/if), 
cos(27r/if), sin(27r/2fX an d cos(27r/2/). Again, we assume 
that the noise decreases quickly with time before merger, at 
t = 0. This model is designed to answer our Question 2 
above. In this case, we find that parameter errors are corre- 
lated only with unique frequency components and the constant 
signal, all the way to tf > O.l/J 1 - The model thus demon- 
strates how components associated with very different vari- 
ation timescales can decouple from each other. Moreover, as 
for the first toy model, we find that marginalized parameter er- 
rors effectively stop improving past a finite time before merger 
(Question 1), which is simply related to their respective fre- 
quencies. As a result, a nonzero residual error remains again, 
even though the signal-to-noise ratio continuously increases 
near merger. 



3. Final toy model 

In our final toy model, we insert a few additional fea- 
tures essential to a realistic LISA data-stream. Firstly, we as- 
sume 5 low-frequency harmonics, 1, sin(27r/i/), cos(27r/if), 
sin(47r/if), (sin47r/if), with unknown amplitudes. We also 
include a high frequency carrier signal with known fre- 
quency, f2 ^> fu but unknown amplitudes in sin(27r/2f) and 
cos(27r/2fX for a total of seven free parameters. Secondly, 
we note that the LISA system is equivalent to two orthogonal 
arm interferometers with both detectors measuring polariza- 
tion phases simultaneously (which correspond to the real and 
imaginary parts of the amplitude modulation, § IV). There- 
fore, the signal is comprised of 4 simultaneous data-streams. 
We incorporate this feature by assuming 4 measurements (i.e. 
4 corresponding Fisher matrices) of the signal with 4 given 
phase shifts (ip* 1 , tp? , ip° 2 , ipf ; 1 < ;' < 4) so that 

hit) = cq + s\ sin(27r/if + Lp s i , ) + c l cos(2Trfit + (p° i 1 ) 

+52 sin(27T fit + ipf) + C2 COS(27T fit + iff) 

+iio sin(27r/ 2 + ci cos(27r/ 2 0, (83) 

In this case, we find that 4 principal components improve 
quickly at late times. As in our second toy model, the high 
frequency parameters decouple from the slow frequency ones, 
except at very late times when ff > O.l/j 

This final toy model allows us to answers all of Questions 1- 
4 as follows. 

• Answer 1: Four out of 5 slow principal components 
of the error ellipsoid are quickly improving with time, 
while one of them stops improving at tf < 0.1/j" 1 . 
Therefore, any parameter with a large projection along 
this one poor principal component will stop improving, 
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while parameters nearly orthogonal to it will keep im- 
proving quickly. Thus, both statements made in relation 
to Question 1 above can in fact be correct, depending on 
the connection between a given parameter and the poor 
principal component. Typically, we expect marginal- 
ized parameter uncertainties to evolve as (S/N)' 1 for 
ff > 0. 1 j\ 1 . For smaller ff values, closer to merger, 
they would continue to improve, albeit with a shallower 
slope. 

• Answer 2: We find that the introduction of additional 
high frequency components does not change the evo- 
lution of original parameter estimations as long as the 
time-to-merger is larger than a fraction of the time pe- 
riod of the additional high frequency components. 

• Answer 3: As the signal-to-noise ratio increases quickly 
at late times, rapidly evolving parameter error combi- 
nations are given by the principal components of the 
error ellipsoid corresponding to the final situation at 
merger. With 4 data-streams, there are 4 such best prin- 
cipal components. Analogously, for the LISA ampli- 
tude modulation given by eq. (25), we expect that the 
2 polarization phases for the 2 beam patterns at ISCO 
can be best determined: (1 + cos 2 6*ail)-F+' /7 (^isco) and 
cos9 NL F^"(Q,isco)- (In terms of ecliptic angular vari- 
ables, these are the real and imaginary parts of the com- 
bination given by eq. (34).) 

• Answer 4: The widths of error distributions for slow 
parameters do not change significantly as long as ff > 
0.1/f 1 . During this final stretch of time before merger, 
however, one of the principal components stops improv- 
ing and the major axis of the error ellipsoid freezes. 
Since the physical parameters can be considered to 
be randomly oriented with respect to the ellipsoid 
axes, distributions of marginalized errors suddenly start 
broadening for ff < 0.1/f 1 , with a worst case relative 
orientation leading to very little improvement and a best 
case relative orientation corresponding to a scaling with 
(S/NT 1 . 



B. Implications for LISA 

These simple toy models offer a general interpretation of 
the time dependence of LISA'S parameter estimation errors 
for source localization. The LISA data stream is described by 
N p i = 5 physical parameters, p s i ow , which are not the harmonic 
coefficients themselves but determine these coefficients, gj 
(or conversely, the mode expansion coefficients gj determine 
the physical parameters p s i ow ; see § III). Neglecting Doppler 
phase and spin precession effects, 2/ max +1=9 modes deter- 
mine the signal by eqs. (41,42). In principle, any N p \ = 5 of 
the gj mode amplitudes uniquely determine the physical pa- 
rameters, psiow However, in the presence of noise, each of 
these modes are uncertain and the combination of all modes 
helps in reducing the estimation errors of the p s i ow variables. 



The key implication of our toy models for LISA is that 
the estimation of low frequency gj modes with low \j\ are 
effectively decoupled from the high frequency signal, unless 
the merger is within ~ 0.1 times the cycle time of the fast- 
oscillating signal. We have shown that the HMD of the or- 
bital modulation consists purely of low-order harmonics, with 
|j I < 4. In comparison, the high frequency GW phase has a 
much higher frequency, corresponding to j > 1000, and this 
high frequency signal's cycle time is greater than the time 
to merger throughout f > fisco Hence, physical parameters 
Psiow will remain decoupled from parameters pf ast , all the way 
to ISCO. This finding is independent of details of the wave- 
form and the modulation, in agreement with the results of 
Ref. [13] which show that decoupling occurs independently 
of the details of the h c (t) signals, including the modified in- 
spiral waveforms of alternative theories of gravity. In terms of 
post-Newtonian expansions, only terms above second order 
have cycle times as large as the cycle time of the amplitude 
modulation. These terms are responsible for the small cross- 
correlations of the two sets of parameters found by Ref. [7]. 

We have not considered spin precession effects, but Vec- 
chio [17] and Lang & Hughes [22] find that spin precession 
effects can help improve the final localization errors by a fac- 
tor of ~ 3. Spin precession cycle times decrease continuously, 
become of order a few days or less during the last week prior 
to merger, and of order hours during the last day of inspiral. 
Therefore, according to our simple models, we expect spin 
precession effects to improve the source parameter estimation 
errors especially during the final two weeks before ISCO. Dur- 
ing that period of time, in the absence of spin effects, param- 
eter uncertainties (especially the sky position major axis and 
the luminosity distance) cease to improve when using only the 
amplitude modulation. 

The best-determined parameters at ISCO are, approxi- 
mately, the independent detector outputs at ISCO, i.e. the real 
and imaginary parts of h{ (p\): d : [ 1 (l+cos 2 9 NL )F+"(S}) and 
d^ 1 cos 0nlF x ' (f2) (see Appendix A 4). These are the 4 in- 
dependent combinations of 5 physical parameters p\ which 
correspond to the eigenvectors of the error covariance matrix 
following the steep evolution cx (S/N)~ l all the way to ISCO. 
We refer to the fifth independent combination, which is or- 
thogonal to these best eigenvectors, the "worst" eigenvector, 
since for this combination, the evolution ceases to improve 
as (S/N)' 1 within ~0.1 x (amplitude modulation cycle time) 
of merger. It is straightforward to obtain this worst combi- 
nation explicitly by using the 4 other eigenvectors and Gram- 
Schmidt orthogonalization (but we have not done this in prac- 
tice). Since the highest frequency harmonic of the slow modu- 
lation is for j = 4, the corresponding cycle time is yr/4. Thus, 
we expect errors will stop improving roughly 1-2 weeks prior 
to merger. Distributions of errors will quickly broaden during 
these final stages of observation before ISCO. Simply scal- 
ing errors with (S/N)' 1 , as in the angle-averaged formalism 
(e.g. [13, 15]), is acceptable if one studies the evolution of 
parameter errors at ff > 2 weeks, or if one only focuses on the 
best case parameter combinations. In general, the exponent 
in the (S/N) scaling decreases as one approaches merger time 
depending on how close the particular combination of angles 
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considered is to the worst combination. 

Our findings for the eccentricity evolution of LISA'S sky lo- 
calization error ellipsoid can also be understood with the sim- 
ple toy models. In fact, we found this behavior to be expected 
for any model signal with relative instantaneous signal ampli- 
tude increasing quickly with time, e.g. t~ a , a > 2. In this case, 
the principal axes of the general parameter error ellipsoid sep- 
arate near tf = 0. There are a limited number of principal errors 
which rapidly decrease to zero near tf = 0, while others "freeze 
out" at a time related to a fraction of the cycle time of the par- 
ticular waveform (<~ 0.ir cyc i e if t\ > r cyc i e ). For LISA, there 
are 5 variable parameters, p s i ow = (dL,0N,4>N,9L,4>L), and es- 
timation uncertainties of 4 combinations of these parameters, 
c/£ 1 (l +cos 2 QNL)F+ H (£l),d.L cosOnlF^ 1 ' (O), improve quickly 
with (S/N)~ l . These combinations correspond to the best 4 
principal axes of the 5-dimensional error ellipsoid. The re- 
maining 5 th principal axis does not improve as (S/N)~\ but 
rather stops improving at a fraction of the last modulation cy- 
cle time. The two dimensional sky position error ellipsoid is 
the projection of the general 5-dimensional error ellipsoid on 
the (9n,4>n) plane. This plane will generally not be aligned 
with the principal axes of the 5-dimensional ellipsoid. In a 
typical case, therefore, there will be a nonzero projection on 
the worst principal component and the sky position ellipsoid 
will stop shrinking along the worst principal component. This 
explains why the major axis, 2a, ceases to improve and the 
eccentricity increases close to merger. 

According to this argument, it is somewhat surprising to 
find that the minor axis, 2b, can stop improving much before 
ISCO. Figure 2 shows that this happens in the worst 10% of 
all cases for randomly chosen source angular parameters. The 
reason for this is that, in some cases, not all rapidly improv- 
ing "best" principal components have a small absolute error 
at ISCO. For example, consider an edge-on binary inspiral 
(cos9 NL w 0). Since two of the quickly improving parame- 
ters are simply proportional to cos 9 NL , the errors will be very 
large for these parameters. Thus, depending on the relative 
orientation of the detector and the source at ISCO, there can 
be large absolute errors in some cases even for the best combi- 
nations of parameters. In short, both axes of the sky position 
error ellipsoid can stop improving at late times in those cases 
when LISA is oriented in its least favorable direction at ISCO. 



VIII. CONCLUSIONS 

We have developed a new harmonic mode decomposition 
(HMD) method to study the feasibility of using LISA inspi- 
ral signals to locate coalescing SMBH binaries in the sky, 
as the mergers proceed. According to our extensive HMD 
survey of potential LISA sources, it will be possible to trig- 
ger large field-of-view searches for prompt electromagnetic 
counterparts during the final stages of inspiral and coales- 
cence. Our results indicate, for instance, that for a typical 
Z <~ 1 merger event with total mass M <~ 10 5 - 10 7 M Q , a 10- 
day advance notice will be available to localize the source to 
within a lOdeg 2 region of the sky. The advance notice to lo- 
calize the source to a 10 times smaller area of ldeg is < 1 



day for the typical event, suggesting that a wide-field instru- 
ment of the LSST class, with a lOdeg 2 field-of-view, may of- 
fer significant advantages over a smaller, 1 deg 2 field-of-view 
instrument for observational efforts to catch prompt electro- 
magnetic counterparts to SMBH binary inspirals. 

The robust identification of such electromagnetic coun- 
terparts would have multiple applications, from an alter- 
native method to measure cosmological parameters to pre- 
cise measurements of merger geometries in relation to host 
galaxy properties [8, 15]. If such electromagnetic counter- 
part searches can be implemented effectively and successfully, 
LISA could become an extremely valuable instrument for as- 
trophysics and cosmology, beyond the original general rela- 
tivistic measurement goals. Given the advance warning time 
capabilities established here, effective strategies for electro- 
magnetic counterpart searches, including the concept of par- 
tially dedicating a > 10 deg 2 field-of-view fast survey instru- 
ment of the LSST class, are considered in detail in a separate 
investigation [18]. 
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APPENDIX A: SIMPLE TOY MODELS 

1. Single Frequency Model 

First, let us consider the following simple model with two 
unknowns, cq and c\ , 

/z(f) = co + cicos(27r/®f), (Al) 

where /® = yr" 1 is fixed and assumed to be known prior to the 
observation. We call t the "look-back time" before merger. 
Let us assume that the relative noise continuously decreases 
during the observation and that the differential squared signal- 
to-noise ratio (without modulation) is given by cr~ 2 (t) = t~ 2 
in eq. (50). Here t = is a proxy for the "merger". Close 
to merger, the signal-to-noise ratio accumulates very rapidly. 
We assume that h{t ) is measured in the time interval t{>t> tf, 
where f; is the start of observation, tf is the end of observation 
(i.e. x = Wrger-f, *min = tf, and x max = fi in eq. [50]). We fix 
f; and examine the dependence of parameter estimation errors 
as a function of tf, assuming tf <C t\. 
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Note that, for the signal (Al), the fiducial values (co,c\) 
drop out when calculating the RMS parameter errors Acq and 
Ac; using eq. (50). More generally, this is true for any sig- 
nal which is a linear combination of the unknown parameters. 
All our toy models will have this property and the results pre- 
sented in this section will be general in that respect. 

First, let us substitute (Al) in (49) and (50), and evaluate the 
expected covariance matrix numerically. Figure 8 displays the 
time dependence of marginalized parameter errors and princi- 
pal errors. The plots show that the parameter errors all de- 
crease with the signal to noise ratio when the look-back time 
before merger is large. However if the end of the observation 
is within a certain critical time to merger, tf < t c , only one 
principal component follows the signal-to-noise ratio. Fig- 
ure 8 shows that f c ~ 0.1 yr. The start of the observation in 
Figure 8 was fixed at t\ = 5 yr. 

It is also interesting to examine what happens for general 
total observation times, do errors stop improving within some 
time f c before merger? If yes, how does f c depend on the two 
timescales t\ and ? We examine this question numerically, 
substituting (Al) in (49) and (50) and now varying both tf/f^ 
and h/fZ l . Let us define the critical end-of-observation, f c , 
as the time when the marginalized squared parameter error 
is first within a factor of 2 of its final value. Figure 9 plots 
the result for the two parameters. Figure 9 shows that t c is 
determined by fq} for large t l7 but becomes f;-dependent fol- 
lower ti values. In the limit t{ <C f^f, the critical look-back 
time is independent of f^ 1 , it becomes a constant fraction of 

Note that, in the limit of an observation extending up to 
merger, at t = 0, the signal becomes h(0) = cq + ci and it has 
infinite instantaneous signal-to-noise ratio. Therefore, this is 
the best combination of parameters for which the scaling of 
errors can follow (S/N)' 1 all the way to t = 0. The worst com- 
bination is cq — c\, which stops improving before t = 0. 

For this simple model, the origin of these features can be 
understood by analyzing the principal errors and the marginal- 
ized errors in the error covariance matrix. For this purpose, we 
present an analytical algebraic solution to this problem. To 
simplify the equations, let us set the time-scale to fZ l /(2tt). 
In this case the Fisher matrix (50) is 
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The integrals can be evaluated analytically, 
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FIG. 8: Marginalized parameter errors (top) and principal errors (bot- 
tom) for the single frequency model. The green curve shows the scal- 
ing with inverse squared signal-to-noise ratio, (S/N)~ 2 , for reference 
on both plots. A total observation of ti = 5 yr is assumed. Marginal- 
ized errors follow the signal-to-noise ratio for large tf, but they stop 
improving within tf < t c ~ 0. 1 yr from merger. Only one eigenvalue 
scales with the signal-to-noise ratio near merger. 



a. Long Observations (/ ffi <C fi) 

Here, we assume that the signal has been measured for 
a very long total time and we concentrate on the effects of 
changing the end of the observation time, tf, near merger. 
Therefore, we take the limit f; — > oo, for which 



r, 7 (Mi) = 



2f£l + Si(f) 



cos(t) 

cos(2f)+l 
It 



+si(0 

+ Si(2f) 



(A3) 



where Si(x) = ^^dx is the sine integral. 

In the next two subsections, we find the limiting behavior 
of marginalized and principal parameter errors in two different 
limits: <C h and f; <C fj* t respectively. 



tt/2 
tt/2 tt/2 



^ + Sife) 

cos^)+l +si(2ff) 



(A4) 



We consider the case of a total observation time which is not 
negligible compared to a cycle time, f^ 1 , i.e. tf <C h. We next 
examine two possible cases, fZ l <C t f <C 1 1 and tf <C /X 1 <C h, 
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errors (i.e. diagonal elements) of (A6) become 
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FIG. 9: Critical look-back time, f c , at which parameter errors stop 
improving. Here f c is defined as the time at which marginalized 
squared errors are within a factor of 2 of their final values for the 
first time. 



separately. 

First let us assume that the merger is still far away in time 
in units of a cycle period (/© <C ff <C h). We substitute (A4) 
in (49) and expand r~'(ff) into a ff 1 series: 



sin(2f r ) . cos(2r r )-l 

~^5T + — S — 




(A5) 



Equation (A5) gives the large ff behavior of marginalized er- 
rors and correlations, which can be compared to Figure 8 in 
the appropriate regime, ff > 1 yr. In this case, to leading or- 
der, all of the squared errors scale with ff, which is the scaling 
of the inverse squared signal-to-noise ratio, (S/N)~ 2 , for our 
noise model. 

Next, let us examine the case when the end-of-observation 
time is close to merger, i.e. ff <C /© <C f;. Now, taking the 
inverse of the matrix and expanding into a ff series around 
ff = gives 



1 + : 



7T 



2tt 2 



1 + fff 



(A6) 



which gives the short timescale behavior of marginalized er- 
rors and correlations. The eigenvalues of T" 1 define the 
squared length of the individual principal axes of the parame- 
ter error ellipsoid, in this case 
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It , 37T f 2 

2 + 16 f f 

7T 2 V 16 7T / I 



(A7) 



Note that, in eqs. (A2)-(A7), time is measured in units of 
fZ l /(2ir). In full units, the squared marginalized parameter 
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For the eigenvalues (A7), we get 



(Sv 2 ) 
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Jj/^ 1 ss 0.1 yr, and it increases 
Equation (A9) shows that one 



Equation (A8) implies that the evolution of the marginalized 
squared error on cq is very flat for small ff, when the second 

term is negligible, i.e. ff -C \J j^rf^ 1 = 0.267 yr, then rises 

steeply (oc f f 3 ). The marginalized squared c\ error is also con 
stant near merger, for ff <C 
oc f f oc (S/N)~ 2 for larger f f 
of the principal errors has a very different time-evolution: it 
has no constant term proportional to t®. Therefore the semi- 
minor axis of the error ellipsoid can decrease continuously 
with the signal to noise ratio. On the other hand, the semi- 
major axis becomes constant for ff -C ^r/© 1 = 0.4 yr. Since 
the marginalized errors are nontrivial linear combinations of 
the principal errors, the constant principal error carries over to 
both marginalized errors and dominates their evolution. All 
of these findings are in excellent agreement with the numer- 
ical results shown in Fig. 8 for ff <C 1 yr and in Fig. 10 for 

h/fg>i. 

It is worth emphasizing that, even if the total observation 
time had been infinite, t\ — > oo, the parameters could not have 
been estimated to infinite precision in this model. It is not 
very surprising if one recalls that in this model we defined 
errors to be infinitely large at infinitely early times (a 2 (t ) oc f 2 ). 
For stationary noise, the contribution of the last cycle to the 
resultant RMS estimation error for a total observation of N cyc 
cycles is 1 / y/N cyc . In contrast, rather than the total number 
of cycles, the typical error during the last cycle dominates the 
determination of noise, for the particular noise model used 
here. 

The main conclusion from this toy model analysis is that 
errors stop improving close to merger, at t c ~ O.lf^ 1 . It can 
be extended to more general noise models, with cr" 2 (f) = t~ a 
and Repeating the calculations for larger a values, we 

find that parameter estimation errors become more and more 
insensitive to very early times, ff <C t ~ t\, and that marginal- 
ized parameter estimation errors cease to improve at some f c , 
which is now an a-dependent fraction of a single cycle time 
before merger. For a > 2, we find that errors increase more 
abruptly at ff > f c , which is consistent with the signal-to-noise 
ratio being a steeper function of time. On the other hand, fol- 
lower a values, parameter estimation errors become more and 
more sensitive to very early times, ff -C t ~ t\. In this case, 
the marginalized parameter estimation errors are again very 
slowly changing for ~ ff < t c , but the approximate time t c 
at which parameter errors stop decreasing will be primarily 
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determined by fi, rather than by the cycle period . The 
transition at ff > f c is not as abrupt, but extends to several cy- 
cles. The a = case corresponds to a stationary instantaneous 
signal-to-noise ratio, with errors scaling slowly as 1 /\A;-ff. 
This case is irrelevant to LISA inspiral signals, which have 
a ~ 2 to a good approximation for lday < t < t\ in the rele- 
vant range of SMBH masses. 



b. Short observations (t\ <C ) 

Let us now examine the opposite limiting case, where the 
start of observation time is already within the final cycle be- 
fore merger. This is relevant to LISA signals, since the ob- 
servation time of SMBH inspirals is often below a full year, 
especially for (1 +z)M > 4 x 10 6 M o . 

We again restrict ourselves to the case with a total observa- 
tion time that is non-negligible, i.e. ff -C f;. Using time units 
of /m 1 /(27r), expanding (A3) into a series of both f; and tf/tu 
we get 




(A10) 



Equation (A10) gives the parameter estimation covariance 
during the final stages of observation before merger for small 
total observation times. In this case, the final errors strongly 
depend on the total observation time. The errors reach their 
final values when the second term becomes negligible in 
eq. (A10). To leading order, this happens at t c ~ ti/3 for 
both parameters, independently of the cycle time, . Equa- 
tion (A10) approximates well the fj dependence of f c shown in 
Fig. 9 for t-Jf^ < 0.2 



2. Double Frequency Model 

Now consider a more elaborate model with five unknowns 
c , si, ci, sio, and c w : 
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FIG. 10: Marginalized parameter errors (top) and principal errors 
(bottom) for the double frequency model. The green curve shows the 
scaling with (S/N)~ 2 for reference on both plots. A total observation 
time ti = 5 yr is assumed. Marginalized errors follow the signal-to- 
noise ratio for large tf values, but they stop improving after ts < 0. 1/, 
for both frequencies. By comparing the two plots, it is clear that high 
frequency component errors decouple and that they are determined 
by two corresponding eigenvalues in the bottom panel. 



h{t) = co + si sin(27r/if) + cicos(27r/if) 

+sio sin(27r/ 2 f ) + ei cos(27r/ 2 f). (All) 

Here, the signal is comprised of two different characteristic 
frequencies, f\ and fa, for which we assume f\ <C f%. More- 
over we assume that f\ and fi are fixed and known prior to the 
measurement, e.g. we take f\ = lyr" 1 and fa = 10yr -1 . We 
again assume an observation in the look-back time interval 
h > t > t{ and take the average instantaneous signal-to-noise 
ratio to increase as a(t)~ 2 = t~ 2 . 

Let us substitute in (49) and (50), and evaluate the expected 
covariance matrix numerically. Figure 10 displays the results. 
As in the previous model, these plots show that all parame- 
ter errors decrease with signal to noise ratio until the last cy- 
cle and all marginalized errors stop improving beyond some 



nonzero residual error at late times. Thus, the general trends 
shown in Fig. 10 are very much similar to the ones in the pre- 
vious simple model (Fig. 8). Again, contrary to the standard 
1 / ^/A'cyc expectation, the error during the last cycle domi- 
nates the total error of the accumulated signal. Moreover, 
comparing Figs. 8 and 10 shows that the presence of addi- 
tional independent high frequency degrees of freedom practi- 
cally does not modify the evolution of marginalized parameter 
errors associated with low frequency components, if f; > /f 1 . 
During the final cycle, the error ellipsoid becomes "thin" and 
the narrow dimension will not be aligned with any of the pa- 
rameters. As a result, this bad principal error dominates each 
of the marginalized parameter errors at late times. (Note that 
the start-of-observation time in Figure 1 is f i = 5 yr.) 
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6h(t) = Sc + 5si sin(27r/ 1 t) + 8c x cos(27r f x t) 
+Ss 2 sin(27r/ 2 f) + Sc 2 cos(2tt f 2 t) 
h = 10/, 



0.1 r 



0.01 



0.001 




<<Hj) 
(fc?) 
<&?) 

<5d) 



0.01 



0.1 



10 



<fti(t) = Sc a + 5s 1 sm(2nfit) + 8c x cos(27r/^) 



0.1 



0.01 



0.001 




0.01 



FIG. 11: Critical look-back time, r c (as in Fig. 9), at which marginal- 
ized parameter errors stop improving. Top: Only (co,si,ci) are al- 
lowed to vary, using the prior (jiojCio) = (0,0). Bottom: All 5 pa- 
rameters (co,ii,c'i,.sio,cio) are determined from the observation. For 
ti > /f 1 , estimations of low frequency parameters (co , Si , ci) stop im- 
proving at t c ~ 0.1/f 1 , while improvement for high frequency pa- 
rameters occurs all the way to t c ~ 0. l/^ 1 ■ 



The critical look-back time, f c , at which this happens is dif- 
ferent for the different frequency components. The top panel 
in Fig. 10 shows that t oi ~ 0. 1/; approximately for both sets of 
components (si,c\) and (sio,cio), where /, denotes the corre- 
sponding frequencies /i = 1 yr" 1 and = lOyr" 1 , respectively. 
The bottom panel in Fig. 10 shows that the principal errors 
separate in three groups. There is one best eigenvector that 
improves continuously until the end, two that stop improving 
near f cl ~ 0.1/i and two that stop improving at f c2 ~ 0.1/2- 
The high frequency parameters (sio, cio) totally decouple from 
the two worst principal components, (vo,Vi), and, as a result, 
decouple from the low frequency parameters (co , s \ , c i ) which 
are primarily determined by (vo, vi). 

As for our previous model in § A 1, the critical look-back 
time is generally different for different t\ values. The bottom 
panel in Figure 1 1 shows the time f c at which the squared er- 



rors first double, as a function of fi/Zf 1 , as in Fig. 9. Fig. 11 
justifies the rule-of-thumb scaling t ci ~ 0. \f\ if the observation 
time is at least one cycle period, /f 1 . 

The central question for the present analysis is how sensi- 
tive is the time evolution of low frequency modulation errors 
to the presence of high frequency components. We can ex- 
amine this question by computing the critical look-back time, 
t c , when the high frequency terms are totally neglected. The 
top panel in Fig. 9 shows that, if one limits the parameters to 
(co,si,ci), and the total observation time is not smaller than 
the long-period cycle time, ~ /f 1 , the resulting f c value for pa- 
rameters co and ci is unchanged at the few percent level. How- 
ever, if the high frequency components are introduced, the S[ 
error evolves differently since it asymptotes already at much 
larger ? c values (~ 0.1/f rather than ~ 0.03/_i). The reason 
is that, for small t, with a noise level decreasing quickly, the 
corresponding function s\ sin(27r/i?) « 2wf\Sit is linearly in- 
dependent of, and thus uncorrelated with, the functions co and 
C[ cos(27r/if ) which are both constant to first order. Hence, if 
there are no more unknowns than (co,*i,Ci), then cq and C[ 
are correlated while S[ is decoupled and can be determined 
independently of the other parameters. However, if we add 
any parameters which are not constant for t <C /f 1 , then s\ 
becomes correlated with those. This is exactly what happens 
in the bottom panel of Fig. 9, when considering the high fre- 
quency modulations: the estimation on s\ becomes limited for 
t < tci ~ 0.1/f 1 due to the correlations with and cio Quite 
similarly, if one introduces any other low-frequency function 
that is not constant to first order, like 52 sin(47r/if), then the 
correlations with this parameter will limit the improvement of 
estimation errors for s\ at t c ~ 0.1/i, even when neglecting 
the high frequency components. As we shall see, this is the 
case for LISA: there are generally more than one sin and cos 
low-frequency modes. In this case, the evolution of estima- 
tion errors for low frequency parameters can be obtained with 
the high frequency modes (like siq and cjo) priored out. This 
justifies our simple intuition: once the signal is decomposed 
into different time-scale components, the parameter estima- 
tion problem becomes separable and the evolution of param- 
eter errors corresponding to different such time-scales can be 
estimated independently from each other. 

Rather than going through an analytical derivation as in 
§ A 1 , we answer one remaining question here: what combina- 
tion of the original parameters (co,*i,ci,sfio,cio) corresponds 
to the best principal component, Vo, which can be determined 
extremely accurately at late times, f f — * 0? At t = 0, the noise 
drops to zero. Therefore, the quantity we can measure us- 
ing the t = information is simply h(t = 0). Looking back at 
eq. (All), this is co + ci +C2- It will be interesting to look for 
similar "best determined combinations" of physical parame- 
ters for the case of the LISA'S realistic signals. 



3. Four data-stream models 

For our final toy model, we insert additional features of a 
realistic LISA data-stream. We consider five low frequency 
unknowns, cq, si, C\,S2, Q, and a high frequency carrier signal 
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with additional unknowns sio, and ck>. Moreover we consider 
the simultaneous measurement of four data-streams. The sig- 
nal is 

h(t) = co+S\sm(2Kf\t + t^) + CiCOz{2'Kfit+(pf) 
+S2 sin(27r/i t + <pf ) + c% cos(27r/i t + tpf) 
+iiosin(27r/2f) + ciocos(27r/ 2 f), (A12) 

where ^ l Sl ' C2 ' 12 (/ = 1 . . .4) are fixed at a priori randomly cho- 
sen numbers defining the relative phases of the various modes 
which are being simultaneously measured. We compute in- 
dependent Fisher matrices for each four set of ^ 1,,Sl > C2 ' i2 . We 
assume that f\ <C fi and that f\ and / 2 are fixed and known 
prior to the measurement. We choose / 2 = 10/i and find the 
evolution of marginalized errors and principal errors in two 
limits: 
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(i) neglecting cross-correlations with the high frequency 
parameters by assuming a prior Ssw = Sew = 0, and 

(ii) accounting for these high frequency parameters. 

We again assume an observation in the look-back time inter- 
val ti > t > tf and take the average instantaneous signal-to- 
noise ratio to increase as er(f)~ 2 = f" 2 . 

The results for these models are shown in Figure 12. Th 
marginalized errors (top) and principal errors (bottom) are 
shown for both cases (i) and (ii) above. The figures show 
that, in agreement with our previous model, uncertainties on 
the low frequency parameters are not affected by the high fre- 
quency parameters, except during the final 0.1 cycle time of 
the high frequency component, O.l/J 1 . The figures also show 
that the four principal components of the error ellipsoid im- 
prove quickly at late times. 

Marginalized parameter errors improve quickly if they have 
negligible projection on the bad directions of the error ellip- 
soid. As a result, our expectation is that errors will typically 
not stop improving abruptly, but that there will be a shallower 
evolution in the final two weeks. In the worst case for a given 
parameter, if it is aligned with the bad ellipsoid principal com- 
ponent, it will stop improving near merger. In the best case, 
if the parameter is orthogonal to the bad ellipsoid principal 
component, it will improve quickly throughout the final days 
of inspiral. Therefore, we understand that the distribution of 
errors broadens for tf <C 0. lf[ . 
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FIG. 12: Marginalized parameter errors (top) and principal errors 
(bottom) for the four data-stream model. Pairs of curves with the 
same line style show results for cases with five and seven param- 
eters. The extra two parameters correspond to high frequency (f£) 
components, which affect errors on the other parameters through cor- 
relations only slightly (factor of < 2) if tf < O.iyj . The green curve 

shows the scaling with inverse squared signal-to-noise ratio, (S/N)" 2 , 
for reference on both plots. A total observation time t\ = 2yr is as- 
sumed. Marginalized errors follow the signal-to-noise ratio for large 
tf values. Four principal errors scale with the signal-to-noise ratio 
near merger. 



4. Best Determined Parameters 

In the previous section, we have shown that, if the noise 
decreases quickly like f 2 near merger (at t = 0), the best- 
determined parameters are the eigenvectors of the error co- 
variance matrix that improve with (S/N)~ l . Near merger, 
these are the independent detector outputs at t = 0. In the 
case of LISA inspirals, the observation only extends down 
to ISCO. In this case, the best determined combination of 
physical parameters p\ at ISCO are the real and imaginary 
parts of h{ (p{). To prove this, we have to show that these 
are uncorrected and decrease with (S/N)~ l . The functions 



/i'(f) and h (t) are uncorrected by construction, since they 
correspond to the two independent Michelson detector out- 
puts (see § II B and Cutler [19]). The real and imaginary 
parts of one of the detectors, dih\(t) and 37i', are uncorre- 
cted since they are the coefficients of the high frequency car- 
rier, sm<j)cw and cos (f>cw, for which correlation over one <j>cw 
cycle (during which the detector noise is approximately con- 
stant) is zero. Another way to see this is to focus on the 
real part in the definition of the Fisher matrix (57), which is 

expressed as the integral of ^t\d a h{ {t)di>h{ (£)]• The term 
in brackets is purely imaginary for the cross correlation of 
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$th\(t) and ^sh\, hence the real part is always zero. There- 
fore, the correlation matrix for 3lh{ ,$sh{ is diagonal. For 
diagonal terms, the derivatives are 1 and the integrals be- 
come simply J a~ 2 dt, which is exactly (S/N) 2 . The RMS es- 
timation uncertainty of 3lh{ (j>\) and $sh{ (p\) follows the 
(S/N) -1 all the way down to ISCO. These best combinations 
are d^(l + cos 2 9 NL )Fl'"(Cl) and cos 9 NL Flf'(Sl). 

The evolution of an arbitrary combination of angles will be 
determined by the projection of this combination on the co- 
variance matrix eigenvectors. A linear combination of good 
eigenvectors leads to similarly quick improvement of errors 
with (S/N) . However, as soon as there is a nonzero projec- 
tion on the fifth eigenvector, the estimation uncertainty will 
stop improving at <~ 0. ir cyc i e which, for the highest j = 4 har- 
monic, is between 1-2 weeks. 



APPENDIX B: ANGULAR VARIABLES 

Here we define the relative angles 9 NL and 4>nl, using the 
polar angles (9 N ,(j) N ) and (Ql,4>l) and the corresponding unit 
vectors N and L. 



Let us write a rotation around z and y as O z (<f>) and O y (9), 
respectively. Then, z = O y (-9 N )O z (-<f)N)N and we define 

sin(fyv L )cos(<^v L ) \ 

sin(9 NL )cos(<j> NL ) = O y (-e N )O z (-4> N )L. (Bl) 
cos(9 NL ) J 

This uniquely defines 6 NL and <f) NL , which correspond to the 
relative latitude and longitude, respectively. More explicitly, 
we get 



we get 

9nl = arccos(N • L) = 



(B2) 



= arccos[sin^Arsin6 | £Cos((/)£-0Ai) + cos6'AfCos6'£] , 
}NL -\ c/> otherwise ' (B3) 



where 

/ cos 9 N sin 

tpo = arccos — — 

\ sinfl NL 



>l cos(cf) L - <j) N ) - sin 9 N cos 9 L \ 

sin0 N L / ' 

(B4) 
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